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ABSTRACT 


No theoretical technique exists yet that enables the engineer to 
completely determine the stability of practical complex spacecraft at- 
titude control systems. In this dissertation, a stability approach is 
developed that utilizes a computer simulation in conjunction with theo- 
retical techniques to locate the potentially troublesome regions of state 
space and explore their acceptability. It is shown in Part I that for 
general physical systems, x = 7(x) , all limit cycles are associated 
with equilibrium points that lie "inside" the limit cycles. Thus, if no 
equilibrium point exists in a region, then no limit cycle is contained in 
that region. Special properties of limit cycles and equilibrium points 
are discussed and illustrated. In Part II, equilibrium points are found 
analytically for spacecraft attitude control systems. A method is devised 
to determine stability near these equilibrium points. Then, instead of 
utilizing a uniform or random n-dimensional grid of computer initial 
conditions, a grid is concentrated around a few calculated points. Sev- 
eral practical examples of actual systems are presented. When control 
systems are viewed in this light, a high degree of confidence is gained 
from far fewer computer runs than previously needed. 
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EQUILIBRIUM ANALYSIS AND ITS APPLICATION 
TO ATTITUDE CONTROL SYSTEMS* 

by 

William N. Redisch 
Goddard Space Flight Center 

1. INTRODUCTION 

No theoretical technique exists yet which enables the engineer to determine the stability of 
practical complex spacecraft attitude control systems. Linearized models and much simplified 
nonlinear approximations of a system are studied on paper, but the resulting analyses are 
nearly always unsatisfactory for describing the stability of the actual system. Hence, the anal- 
ysis of the complete system relies on the use cf detailed computer simulations. 

This approach, necessary as it is, leaves much to be desired. Firstly, haw accurate is the 
simulation? There are oddities in physical hardware which are not known, understood, or 
sufficiently determinable to generate a mathematical model that can be simulated. Other meas- 
urable non-primary characteristics are labeled "higher- order effects" and, to keep the com- 
plexity cf the simulation reasonable, are not simulated. (For example, secondary emission of 
x-rays in a radiation belt which causes spurious noise effects on field-effect transistors.) How- 
ever, the limitations of a computer analysis are acceptable if the engineer takes enough pains 
to ensure that the differences between real and simulated performances are well below the 
error tolerances specified for the actual system. 

A more important shortcoming of computer analysis for complex high dimensional attitude 
control systems is that an infinite number cf stability runs from an infinite number of 
initial conditions obviously cannot be made. Hence, it is possible for a limit cycle, which is 
considered unstable behavior (for the purposes of this dissertation), to exist in state space and 
never be discovered until after the launching cf the spacecraft. That is, the limit cycle-exists 
in a region never entered by the simulated state vector, because every initial state produced a 
state-space trajectory that did not pass through the troublesome region. There are many ex- 
amples of this occurrence, one of which is discussed in Chapter 9. 

In this dissertation, a method is developed that employs a complete computer simulation 
in conjunction with theoretical techniques so as to yield a very high degree of confidence in a 
stability analysis. That is, after a finite number of runs, the simulated state vector will, with 


■'Taken from a dissertation submitted to the Faculty of the Polytechnic Institute of Brooklyn in partial fulfillment of the require- 
ments for the degree of Doctor of Philosophy (Electrical Engineering) 1967. 
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high probability, pass through all troublesome regions in the overall state space of interest. 

Now, if a system is truly globally unstable or stable, then the response either diverges from or 
converges to the designed final state (call it the origin) from anywhere in state space. This 
type of behavior is easily observed on a computer. It is the more nebulous limit cycle or mean- 
dering patterns that are difficult even to define from a stability point of view, but are unwanted 
from an attitude control designer's point of view, that elude the computer analysis. Thus, a 
theory is needed to indicate the locations of all the potentially troublesome regions. The com- 
puter is then used to explore these regions and find out whether they are acceptable. Also, 
once a method of finding trouble spots is developed, insight into a system is gained, and it is 
sometimes feasible to design the system so as to eliminate these troubles. Such an approach 
is illustrated in Chapter 9. 

These ideas are not limited to attitude control systems, but apply to systems in general. 
Hence, this dissertation is divided into two major parts. In Part I, the ideas are put on a math- 
ematical footing applicable to all systems via the linking of system trouble areas to mathemat- 
ically defined limit cycles. It is shown through general functional relationships that the trouble- 
some regions are always associated with equilibrium points. For second-order systems, it is 
well known that singular points play an important role in determining the behavior of solutions. 
Here, it is shown that for general n th -order systems, all limit cycles are associated with equi- 
librium points which lie "inside" the limit cycles. Thus, if no equilibrium point exists in a 
region, then no limit cycle is contained in that region and, engineering-wise, the complete phys- 
ical system does not meander about in this region. The equilibrium points, are found analytically 
and are usually finite in number. Now, instead of a uniform or random n-dimensional grid of 
computer initial conditions, a grid concentrated around the calculated points is used. 

Part llpresents an engineering application of the equilibrium analysis technique for at- 
titude control systems. Necessary equations are developed, and methods for finding, investi- 
gating, and interpreting equilibrium states are discussed with ample illustrations from actual 
systems. When systems are viewed in this light, a high degree of confidence is gained from 
several hundred computer runs, where previously the same level cf confidence could be ob- 
tained only after several hundred thousand runs. 
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Part I 


MATHEMATICAL DEVELOPMENT FOR GENERAL SYSTEMS 

2. GROUND RULES AND DEFINITIONS 

As a first step, some ground rules and definitions are established. Only systems adequately 
represented by deterministic models in aphysical sense are considered. That is, the system equa- 
tions have bounded solutions that areunique withrespect to any initialconditions. The parameters 
within the systemaretime-invariant. Theentire systemisdescribable by means of a finite number of 
finite-order, autonomous, ordinary differential equations, generally non-linear in nature. 

A useful theoretical tool often used in the development is the formulation of the mathemat- 
ical model that represents the system in state-variable form. Thus, if x L ft) is a state 
variable (i = 1, 2, • * * , n) for an n- dimensional system, and t represents continuous, non- 
negative, real time, then the system is described by the following n equations: 


dx x (t) 


dt - f l ( X 1 ( t )> x 2 (t)> " 

• - X n(0) 


dx 2 (t) 



dt " f 2 ( X 1 (<0* x 2 (t)* • • 

• - X „(t>) 

(2.1) 

d^CO 



dt - f „ ( X 1 (t)> X 2 (*)■ •’ 

' - X n (O) - 



where n is a finite positive integer. A shorthand derivative notation used is x = dx/dt, 
x - d 2 x/dt 2 ,etc., or x (i) = d i x/dt i where i = 0, 1, 2, • • and when i = 0 then xW = x . 
Thus, Equations (2.1) are equivalent to 

i i (O = f i ( x i (t). x 2 (t). ••• ,x n (t)} i = (2.2) 

It is also convenient to introduce a vector formalism by defining 



(2.3) 
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and 


hi (*)\ 

' f 2 (*) 

l f »W 


(2.4) 


such that Equation (2.1) can be rewritten as 

^ = T(x(.,) • (2.5) 

Equation (2.5) satisfies the imposed uniqueness constraint, that is, it has a unique con- 
tinuous solution x(t) for every initial condition %(0) = k, ,if 


x - x Q min (a, bM) , 1 1 1 -4 b ; 


( 2 . 6 ) 


||f(x)|| < M ; (2.7) 

and f(x) satisfies a Lipschitz condition in the domain ||x-x 0 ||< a. (Reference 1, Chapter 6, 
and Reference 2 r Chapter 3.) M, a and b are the constants which define the uniqueness 
domain. The function 7(x) satisfies the Lipschitz condition, provided that a constant K is 
found such that 


|f(xi) - f(x 2 )|| < K||'Xj - X 2 | 


( 2 . 8 ) 


where || || represents the norm of a vector. 

Another form in which a system is sometimes represented is 


x< n >(t) = f(x( n " 1 ) (t), x< n " 2 ) (t), ••• , i(t), x(t)) . (2.9) 

Neither Equations (2.1) nor (2.9) contain time explicitly in accordance with the previous ground 
rules. Equation (2.1) is the general state variable form, but form (2.9) often has the advantage 
of facilitating mental visualization of time-trajectories, and can always be changed to form (2.1) 
by letting x. = x< i-1 >, where i = 1, 2, • • * > n • However, in general, it is not possible to com- 
bine n first-order differential equations to one n th -order differential equation (Reference 1, 
Section 2.3 and Reference 2, Section 6.4). Hence, the state-variable formulation (2.1) is more 
general and encompasses form (2.S). 
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For an n-dimensional system, n > 2, described in the state-variable form, a limit cycle is 
a closed, bounded trajectory repeatedly traversed with time in a positive time period r in 
n-dimensional state space. Several authors (References 3 and 4) distinguish periodic solutions 
from limit cycles, calling limit cycles the periodic motion that neighboring solutions tend to as 
t -co. Here, the two are considered as the same. For the proofs in Chapter 3, if a limiting 
type cf limit cycle is the only periodic solution that exists, initial conditions are chosen to lie 
on the limit cycle such that the solutions are exactly periodic from time equals zero. Thus, 
each state variable 


x i (t) “ x i ( t + t) i = 1 , 2 ,'", n ( 2 . 10 ) 

is aperiodic function of time with a common period r. Each has an average value over the 
period defined by 


X 


i avg 



(A) dA 


1 , 2 , 


( 2 . 11 ) 


In n-dimensional state space, the n average values cf the n state variables define a point which 
shall be called the centroid of the limit cycle. For example, in the second-order system 


x 2 - - co 2 Xj + K , 


( 2 . 12 ) 


whose solution is the limit cycle 

(t ) - [ x 2 (0)/&>] s i n w t + Jxj (0) - K/c<; 2 ] cos wt + K/oi 2 

(2.13) 

x 2 (t) “ x 2 (0) cos wt - "[xj (0) - V" 2 ] sinwt , 


the centroid is at 


X lavg = V" 2 X 2 avg ” 0 ( 2 - 14 ) 

in Xj, x 2 state space, and r = 2n/a, This 
limit cycle and its centroid are shown in 
Figure 2.1. 



Obviously the centroid cf a limit cycle 
lies within the limit cycle where "within" is 


Figure 2.1 - Limit-Cycle Centroid. 
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defined as follows: If there is a closed bounded curve in n-dimensional state space, then its 
projection on any of the state variable axes is a finite line segment. A point such as the centroid 
being within the n-dimensional limit cycle means that 

x^ < x iavg - x iM * i = 1, 2, ' ' ‘ , n , (2.15) 

where x im and x M are algebraic quantities representing the end points of the finite line seg- 
ment projection on the x th axis. 

This definition cf "within" does not lend itself to the normal conceptual visualization cf 
"inside" a non-simple flat closed curve. For example, suppose that the plane projection of an 
n th -order limit cycle is a figure-8-type curve as depicted in Figure 2.2 The projection cf the 

centroid on this plane lies anywhere in the 
shaded region. 

An equilibrium point in n-dimensional 
phase space is defined as a point at which 
the solution to the systems' equations is con- 
stant with respect to time. That is, if a set 
of compatible initial conditions exist such 
that the left-hand sides of Equations (2.1) 
remain zero for all time 0 < t < <» , then these 
sets cf initial state variable values, obtained 
Figure 2.2— Area Within a Limit Cycle. from the solutions to 


f 1 ( X 1 > X 2’ 

• X n) " 0 


f 2 ( X l- X 2’ "• 

, x") = 0 

(2.16) 

f n ( X l- X 2’ ••• 

, xl = o , 



are equilibrium points. Hence, equilibrium points are trajectories which are the points them- 
selves remaining stationary in state space for all positive time. The points may be stable or 
unstable equilibriums in the sense that another set of initial conditions, arbitrarily close to an 
equilibrium point, yields a trajectory that tends toward or away from the equilibrium point with 
increasing time. 

In the next chapter it is shown that all solutions to Equations (2.16) are indeed equilibrium 
points, and that when limit cycles exist for an n th -order system (2.1), then there is one or more 
of these equilibrium points within each limit cycle. Hence, a computer simulation of system 
(2.1) can be used much more effectively by investigating the system around these calculated 
equilibrium points. It is noted that the solution of Equations (2.16) is a straight-forward alge- 
braic problem that does not involve the dependent variable time. In fact, this straightforward 
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but tedious task is often automatically performed by the computer simulation itself, whether 
analog, digital, or hybrid. 

Again it is emphasized that this type of procedure, though not 100-percent encompassing, 
is a powerful one for systems in general. The approach is not tied down to specific (attitude 
control) systems until Chapter 6. In Chapter 8 it becomes clear that attitude control system 
equilibrium points do not readily reveal themselves; methodical procedures must be employed 
to find the more elusive ones. 
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3. THEOREMS AND PROOFS 


In this chapter, preliminary theorems are developed following the ground rules and defini- 
tions of Chapter 2 . The theorems are then used to prove that at least one equilibrium point lies 
within a limit cycle when the limit cycle exists in an n th -order system. 

Firstly, if x. (t) and all its derivatives exist for all t, and 

x.(t) = x.(t+T) i = 1 , 2 , ■ . ■ > n , (3.1) 

then it immediately follows that 


x/ k ) (t) = x/ k )(t+r) k = 0, 1,2, ••• (3.2) 

for all t including t = 0. Thus, the relationship in Equation ( 3 . 1 ) adequately defines a limit- 
cycle behavior for the system (2.2)without the necessity cf stating that the derivatives are also 
periodic. 

Since 


x i< k+1) (t) = ^x.< k >(t) 


( 3 . 3 ) 


fork = 0, 1 , 2 , • • *, and i = 1 , 2 , • • • ,n, then 


/•t+T 

x/ k > (t +t) - x/ k > (t) = J x/ k+1 > (k) elk , (3.4) 

assuming x< k+1 ) has no singularities as dictated by the ground rules. Equation ( 3 . 2 ) implies the 
vanishing of the left-hand side of ( 3 . 4 ); thus, 


J x/ k+1 ) (A.) dk 


TX 


(k + 1) 
i avg 


0 . 


( 3 . 5 ) 


Since t > 0, 


X 


(k + 1) 
i avg 


0 for k - 0, 1, 2, • • • . 


(3.6) 
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In other words, 


x (k+i) 


i a v g 



+1 >(\)d\ 



dx ' 


0 


that is, the derivatives of aperiodic function have zero averages. 
It is also clear that 


(3.7) 


x.( k ) (t ) - x (t + t) - xd k ) ((t + r) + r) 

( 3 . 8 ) 

= x3 k > ( ( (t + r) + r) + 

or 


xd k ) (t) - x (k ) (t + jr) 


(3.9) 


for periodic x i (t), where i = 1, 2, • . . ,n; j = 0, 1, 2,. . .,and k = 0, 1, 2,. . 

Another useful fact is that no derivative of a non-constant periodic function is identically 
equal to zero. This is shown by expanding in a Fourier series, 


CO 

x(t) = x( t + t) = x avg + J~|(a n cosont + b n sinwnt) 

n - 1 


CO 



( 3 . 10 ) 


Hence, 


where 


x (k) = x (k) ( t t 7 ) 


CO 

^ ( A n cos ant +B n sinant) 


{ (-l)< k_ 1 > / 2 (ari) k b 
-l ) k/2 (a n) k a n 

{ (-l)< k+i y 2 (m) k a 
(-l) k/2 (on) k b n 


k = 1, 3, 5,- * • • 

k = 2, 4, 6, • - - 

k = 1, 3, 5, . " 

k = 2, 4 , 6, 


( 3 . 11 ) 


( 3 . 12 ) 
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Now, 


n - — I x^ k ) (t) cos or>t dt 

Jo 

(3.13) 

2 f T 

B n = — I ( t ) sinant dt 

Jo 

Thus, if x< k > s 0, then A n = B n = 0 for all n; hence, a n = b n = 0 for alln, implying that x(t) 

= x = constant, which is a contradiction. 

In Chapter 2, equilibrium points of a system are defined as the points in state space cor- 
responding to constant solutions of the system. Hence, when a system yields a unique solution 
for any initial conditions, then, if an equilibrium point does exist, choosing it as the initial state 
of the system ensures that all the derivatives of the system are zero and stay zero. That is, 
for the n-dimensional system, 


d _ 
dt *<*) 




(3.14) 


if 


£(*e) = 0 


(3.15) 


then the initial condition 


x(0) = x e 


(3.16) 


implies that 


x(k) (t) s 0 k = l, 2, • • ■ (3.17) 

for allt > 0, since x(t) - x e is a unique solution for the system (3.14). 

However, there is not necessarily a real x e that satisfies Equation (3.15). For example, 
consider the simple unique system 


_d / x i (t) \ _ / X 2 Ct) 
dt \ X 2 C t ) / " \ 1 


(3.18) 
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which yields the solution 



x i (0) +x 2 (0) t 
x 2 (0) + t 


(3.19) 



as depicted in the state space shown in Fig- 
ure 3.1. 

This system has no equilibrium points, 
since 



has no solution, that is 


(3.20) 


i 2 = 1 / 0 . (3.21) 

But, Equation(3.21) also implies that no limit cycle exists, since by Equation (3.6), x 2avg is zero 
€orlimit-cycle behavior, thus requiring x 2 (t) to be zero at some times , t 2 , • • . , which 
Equation (3.21) does not allow. Thus, the same condition that prevents the existence of an 
equilibrium point also prevents the existence of a limit cycle. 

The previous example is a very special one where f 2 (x lt x 2 ) is actually independent of 
the state variable x. In the next chapter, examples are presented that illustrate systems 
that have limit cycles and yet have no equilibrium points. With each example is an ex- 
planation of those ground rules in Chapter 2 that are being violated. For the systems described 
in Chapter 2 that have one or more limit cycles as solutions, at least one equilibrium point 
exists. This is implied by the proofs that follow shortly, which show that at least one equilib- 
rium point lies within a limit cycle whenever the limit cycle exists as a solution to the system 
equations. 

One additional area of interest is noted here. Instead of working in Euclidean state space, 
it is sometimes more convenient to investigate the topological behavior cf solution trajectories 
in different spaces. For example, if one or more coordinates are cyclic (such as an angle 
determined by only modulo 2n), then a cylindrical or toroidal rather than a planar representa- 
tion may be more suitable. Chapter 8 of Reference 5 is devoted to this subject. The main point 
to note here is that an open curve in a plane can close on itself around the surface of a cylinder, 
thereby acting as a closed limit cycle in cylindrical space, and yet not encircle any equilibrium 
or singular points. Hence, it is re-emphasized that this dissertation deals with Euclidean state 
space only. 
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By definition, an equilibrium point never lies on the trajectory of a non-trivial limit cycle; 
if it did, it would mean that the state cf the system changes with time (a contradiction). Now, it 
is demonstrated how at least one equilibrium point lies within a limit cycle. Instead of dealing 
directly with the general n th -order nonlinear system, the following class cf systems is con- 
sidered in sequence: 

A) Linear systems 

B) Second-order systems, x = f(x, x) 

€) n th -order systems, x (n) = f(x< n_l >, x ( n “ 2 ) ,• . • x ) 

D) General systems, x = f(x). 

A) Linear Systems. 

For linear systems, not only does an equilibrium point lie within a limit cycle, but, also, 
the centroid of the system is an equilibrium point (which obviously lies within a limit cycle). 
The centroid of a system is defined in Chapter 2 as the point x. 

J 1 r iavg 

Consider first the linear system 


n 

(n) (t) = £\-i x (n ' j) (t) 

j=l 


+ C , 


( 3 . 26 ) 


where the a’s and c are arbitrary constants, and 


is a solution. Now, 


x(t) - x(t +t) t > 0 


(n) = 


^ J x (n) (A.) dX 


n 

^ * a n _. x< n -i> (X) +c 
i=l 


dX 



i*t+T /*t+T 

7J x( n ">) (X) dX + t J cdX 


j = i 


( 3 . 27 ) 


( 3 . 28 ) 


n 

E a . x + c . 

n-j • avg 

j=l 
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But, from Equation (3.6), 


*.£? = 0, k = 1, 2, ••• . 

Hence, 

0 - a„ x + c . 

0 avg 

But, from Equation (3.26), the equilibrium value of x is 

0 = a 0 X e + C ' 


Thus, x avg is an equilibrium point. 

For the general state variable formulation of the linear system, 


(3.29) 


(3.30) 


(3.31) 


X i (t) = 


= £■** *■ 


+ a i0 . 


i = 1, 2, • • • , n 


(3.32) 


where the a Lj are arbitrary constants, let x(t) - x(t +r) be a solution. Then, 


/•t+T 

if 

•* t 


x ; - t- | Xj (X.) cE 


_j =i 


+ a . r 


dX. 


But also, by definition, 



aij Xiavg 


+ 


a i0 


= o . 


n 



j=i 


(3.33) 


(3.34) 
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and therefore the x. avg are the x ie , i.e,, the system’s centroid which lies within the limit cycle 
solution is an equilibrium point. 

It is noted that the above proofs hold for any system 

(t) = (x t (t), x 2 (t), • • • , x n (t)) , i ■= 1,2, ---.n (3.35) 


for which 




X 1 avg * X 2avg* 


X 


n avg 


) 


■*r 


f i '(xj (k), X 2 (X), 


< n (M) dk 


(3.36) 


B) Second-Order Systems. 

Consider next the second-order system of the form 

x(t) = f(x(t), x(t)) . (3.37) 

The Van der Pol plus many other well-studied equations are of this class. Assume that for 
some initial conditions Equation (3.37) has a solution x(t) = x(t +r) which corresponds to a 
closed trajectory in the x, x phase plane. The equation of this curve is 

g(x, x,K) = 0 . (3.38) 

Here, as depicted in Figure 3.2, the constant K is determined by the initial conditions that give 
rise to the limit cycle. That is. 


g(x(0), x(0), K) = o . (3-39) 

The slope of the trajectory is expressed as 

dx _ x _ f(x, x) 

sl = dx = l = j ’ (3.40) 

An equilibrium point x = 0, x = x e is given by 

0 = f(0, X e ) (3.41) 


14 


and corresponds to a singular point of indeterminate slope. (In engineering parlance, a singular 
point is a point such that if the system starts there initially, it may mathematically leave that 




X 


Figure 3.2 — Limit-Cycle Trajectory. 


point at various angles. But since the solution 
trajectory is unique, the system remains sta- 
tionary, in equilibrium, at the singular point.) 


The bounded periodic x(t)has minimum 
and maximum values x m andx M , respectively. 
When x( t ) is at one cf its extremums, x( t ) = 0. 
That is. 


g(0, x m ,K) = o 

g(°, x M , R ) = 0 


(3.42) 


Since x avg = 0 lies between the algebraic extremes of x(t ), it is only necessary to show that 


(3.43) 


in order to demonstrate that an equilibrium point lies within the limit cycle. Obviously x e / x m 
and x e / x M since an equilibrium point is not on the limit cycle except for the trivial case when 
the limit cycle is a constant solution. Thus, it lies inside or outside the closed trajectory. 
Along the curve 


x(x, x) = f (i, X) = 0 (3.44) 

in the phase plane, the slope of a solution trajectory is zero, (that is, the solution trajectory 
has a horizontal tangent), as evidenced by Equation (3.40), except (possibly) when x = 0. (The 
locus cf zero slopes x = 0 is used in the Lienard construction method, for example, in Ref- 
erences 5 and 6). The points where the curve 5 = 0 crosses the x-axis, (i.e. “= 0), are equi- 

librium points, since at these points Equations (3.41) and (3.44)become identical. Now, in the 
upper half of the phase plane where x > 0, x(t ) increases with increasing time; hence the solu- 
tion trajectories go horizontally from left to 
right across the zero-slope curve. The op- 
posite is true in the lower half plane, where 
x < 0. See Figure 3.3. Since a limit-cycle 
solution trajectory is abounded, closed, con- 
tinuous curve, it has a horizontal tangent (at 
its crossings of the zero-slope curve) at 
least twice, that is once for x > 0 and once for 
x < 0. The limit cycle therefore encircles at 
least one cf the x e points (where the zero- 
slope curve intersects the x = 0 axis). That 
is, at least one equilibrium point lies inside 
the limit cycle. 
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The above conclusion (that any general second-order system limit cycle encircles a singu- 
lar point) is well known. It can be simply developed via Poincare's theory by showing that at 
least one singularity lies within a simple closed curve since the total Poincark index is +1 (Ref- 
erences 3, 4, and 7). However, the above development lends itself to extension to higher-order 
systems via the following observations. x(t) reaches a maximum value x M at instants of time 
when x = 0 and x < 0, and a minimum value x m when x = 0 and x > 0. The curve x = 0 divides 
the x, x phase plane into two regions, one where x is positive, and the other where x is negative. 
(One of these regions may be nonexistent when a limit cycle does not exist, such as the example 
of Equation (3.18), x = il. However, when a limit cycle solution does exist, then x avg = 0 as 
previously proved. Hence, x is positive and negative at different times, implying the existence 
of phase-plane regions where x is positive and negative.) Furthermore, the x-axis is divided 
into positive and negative x segments, with the boundary points being equilibrium points. Hence, 
at least one of these boundary equilibrium points lies between x m and x H on the x-axis. That is 

x m < X e < x m ■ (3.45) 


C) n th -Order Systems. 

For the class of n th -order systems which yield periodic solutions, 

x<">(t) = f(x("- 1 > (t), x( n - 2 > (t), ••• , x(t>) , ( 3 - 46 ) 

where n >2, similiar reasoning as for the second-order case applies. Firstly, since for 

k = 1 , 2 , • • • f n > 


x ( k ) < x - 0 - x < xd k ^ , 

m e avg M 

it is only necessary to establish that 


(3.47) 


X m < X e < X M (3.48) 

in order to state that an equilibrium point lies within a limit cycle. To establish (3.48) 
for an n th -order system as in (3.46), it is again observed that, even if n is greater than 2, sim- 
ultaneously x = x M when S < 0, and at some other time, x - x m when x > 0. Now, x = 0 is an 
(n - l)-dimensional hypersurface which partitions x, x, 'x, x (4) , • • • , x< n > n-dimensional space 
into two regions; in one, x is positive; in the other, x is negative. Again, both regions exist if 
a limit cycle exists. The equilibrium points of the system are where the x-axis intersects the 
x = 0 hypersurface. Thus, the x-axis is segmented by the x = 0 hypersurface such that equi- 
librium points separate the positive and negative x regions. An example is illustrated in Fig- 
ure 3.4. The projection cf a limit cycle on to the x-axis has its extremums separated by at 
least one x = 0 boundary, and thus at least one equilibrium point. That is, the relationship in 
Equation (3.48) is established. 
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D) General Systems. 

The remainder of this chapter is devoted 
to general systems of the form 

x = T(x) . (3.49) 


Shaded line segment: x >0 

Unshaded line segment: x <0 The approach to the problem of showing that 

at least one equilibrium point lies within a 
Figure 3.4— Positive and Negative ;-Segments. limit cycle solution is via level curves and 

surface gradients, similiar to a geometric 

interpretation of Liapunov functions (for example, in Reference 8, Section 13- For the general 
state variable representation of an n th -order system as in Equation (3.49), the state variables 
x. are not necessarily time derivatives of each other. Hence, in general, 


X . / 

i avg ' 


for any of the i's, i = 1, 2 r • - • ,n. Thus, it is imperative here to prove that 


(3.50) 


< X. < X 


iM 


(3.51) 


for all i , i = 1, 2, • . • s n , rather than the previously simpler task where only one variable, x, 
had to be considered. Here again, expression (3.51) has no equality signs, thus indicating that 
the trivial static case of x. m = x. e = x. M is being ignored. 

First, a general second-order case is discussed since its easy visualization tends to make 
a geometric interpretation significantly clearer. Assume the unique system 


ij (t) = f t (x t (t), x 2 (t>) 

i 2 (t) = f 2 (Xj (t), X 2 (t)) 


(3.52) 


that yields a finite limit cycle solution 

x i (*) = x j (t +t) 


x 2 (t) 

for certain initial conditions x x (0), x 2 (0). 
equivalent system 


(3.53) 

= X 2 ( t+T ) 

If time is eliminated in Equations (3.52), then the 


f 2 ( x i> x z) 

f i ( X l> X 2 j 


(3.54) 
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yields an integral solution which is abounded closed trajectory in the x x , x 2 plane. This limit 
cycle is expressed in parametric form with t as the parameter as in (3.53), or more directly 
by the expression 


g(xj (t), x 2 (t)) = K L c (3.55) 

k l.c. is determined by any point (xj, x 2 ) on the limit cycle, including the initial condi- 
tions (xj (0), x 2 (0)). Since the state of the system is confined to the curve ( 3 . 55 ) for all t, 

g( x , (t), x 2 (t)) - 0 (3.56) 

for the initial conditions x t (0) , x 2 (0). This can be envisioned from a Liapunov point of view 
by taking as a Liapunov function 


v( x i ■ x 2 ) = 

such that for the limit cycle case 

v(x x , x 2 ) = g(x x , x 2 ) = o . ( 3 . 58 ) 

That is, the system trajectory never leaves 
the level curve K L c of the surface V = g. 
See Figure 3 . 5 . 

Two vectors are now introduced: a field 
vector and a gradient vector. The field 
vector is represented by 

& = 

or, for geometric purposes, 

S = Xj x 


g ( x j i x2) > (3.57) 



Figure 3.5— Surface Function g(x lf x 2 ). 



+ x 2 x 2 , ( 3 . 60 ) 


where the x x and x 2 are unit vectors along the positive x x and x 2 axes respectively. Hence, for 
the system (3.52), the field vector 


s ; /f i ( x i (t)- X 2 (O) x i + f 2 ( x i (*)- X 2 (t>) X 2 


(3.61) 
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is geometrically tangent to the solution trajectories including the limit cycle in the , x 2 state 
plane, and points in the direction of increasing time. 

The gradient vector of g(x lt x 2 ) is represented by 


V g = 



( 3 . 62 ) 


or 


Vg 



dg 

<9x 2 x 2 ' 


( 3 . 63 ) 


Hence, the time derivative of g is simply the inner or dot product of Vg and S, That is. 


g - [(Vg, S)] - Vg • s 

( 3 . 64 ) 

- • + JjL ■ 

<9Xj X 1 5 x 2 X 2 



Vg 

to outwards) as the limit-cycle trajectory is 
traversed. Only unique single-valued func- 
tions g(xj, x 2 ) need be considered. That is. 
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one peak (actually an odd number of extre- 
mums)where Vg = 0 that lies within the limit 
cycle. This is shown a bit more rigorously 
by examining g(xj , x 2 ) versus x t for various 
values of x 2 , that is, scanning across a con- 
stant x 2 = K. By Rolle's theorem (Refer- 
ence 10 ), dg/dxj | x ^ =K = 0 occurs at a value 
x j such that x lm 5 a < Xj < b < x 1M , where 

g(a, K) = g(b, K) = K l c . See Figure 3.8. By 
allowing K to vary from x 2m to x 2M it is seen 
that a portion of the curve dg/dx l (xj, x 2 ) = 0 
lies within the limit cycle. Since at the point 
on the limit cycle where x 2 = x 2m or x 2M , 
i 2 = 0 , then 0 = g ^(dg/dx^Xj from Equa- 
tion (3.64). But, except for the trivial case 
of the limit cycle being an equilibrium point, 
Xj / 0 at the x 2 extremum points. Hence, at 
these points, 

= 0 . 

X 2M orx 2 m 
e= K L.c. 




A similiar symmetric argument is pre- Figure 3.8— Section Views of Surface g(x lf x 2 ). 

sented for g-versus-x 2 curves for constant 

Xj values, to show that a portion of the <?g/dx 2 = 0 curve lies within the limit cycle, and that 


<9x 2 


X 1~ X 1M or x lm 
fe =K L.C. 


0 . 
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Now, Vg = o occurs at the points of intersection of the two continuous curves dg/dxj = 0 and 
dg/<?x 2 = 0; at least one such point occurs within the closed limit-cycle trajectory. 

Now the question is to determine if this point where Vg = o, that lies within the limit cycle, 
is indeed an equilibrium point. To do this, a new system is defined as 


. _ ag^.x,) 

x ! ax 2 


. dg( x l' x 2 ) 

x 2 dx x 


(3.65) 


where the function g(xj, x 2 ) is the same solution function to the original system as that given 
by Equation (3.55). Eliminating time from Equation (3.65) yields the equivalent system 


dx 2 dg/dXj 

dXj ~~ <3g/ 3 x 2 


for which the integral solutions are 


(3.66) 


g(xj, x 2 ) - K . (3.67) 

That is, from the definition cf the new system (3.65), all the solutions are closed curves in the 
Xj, x 2 plane, as for example those shown in Figure 3.9. 

Along the limit cycle of the original system and at its equilibrium points, g = 0. Thus, 
from (3.64), 



^ f 2 = o (3.68) 

or 

f 2 dg/dXj 

fj dg/dx 2 ' (3.69) 

Hence, when K - K L c , and at all the equi- 
librium points, the two systems (3.52) and 
(3.65) are equivalent, that is, they yield the 
same unique solutions for the same initial 
conditions. The equilibrium points cf the 
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new system (3.65), by the very definition of that system, are at 


d% _ 
ax 2 


o 


dg _ 
dx t 


0 , 


(3.70) 


or in other words, 


Vi = 0 . (3.71) 

But at least one point for which V g = o lies within the limit cycle g - K L c (for example, the 
point g = K e in Figure 3.9), as previously shown. Since the two systems are equivalent at the 
equilibrium points and on the K L c limit cycle, it is true that for the original system (3.52), at 
least one equilibrium point lies within the K L c limit-cycle trajectory just as it does for the new 
system (3.65) . Thus, the proof that at least one equilibrium point lies within a limit cycle is 
established for second-order systems by this gradient technique. 

These arguments are now extended to higher-order systems. A third-order system is con- 
sidered before finally treating an n th -order system. 

Consider the unique system 

*i (F) = f i ( x i (O- x 2 < t )> X 3 C 4 )) 

x 2 (t) = f 2 ( x ! (t), x 2 (t), x 3 (t)) (3.72) 

x 3 (t) = f 3 (Xj (t), x 2 (t), x 3 (t>) , 

which yields a finite limit cycle 

x i (t) = x £ (t+r) 3 - 1,2, 3 (3.73) 

for certain initial conditions x t (0). In three-dimensional state space, the bounded closed tra- 
jectory is described functionally by the intersection of the two surfaces 


g 1 (x 1 ,x 2 , x 3 ) - K 1LC 

g 2 (Xj , x 2 , x 3 ^ - K 2L c 


(3.74) 


22 



This closed curve has projections in the x 2 , x x x 3 , andx 2 x 3 planes that are described 
functionally by 

h 3 (xj , x 2 j - c 3 

h 2 (x 3 , Xj) - C 2 (3.75) 

W ( X 2 ’ X 3 ) “ C 1 

respectively. The c's and K's are constants. It is also true, for the limit-cycle initial conditions, 
that 


9 


= 




= 0 . 


Introducing the field vector 



and the gradient vectors 



then 


h. = [(vtv.s)] = 0, 


where [( )] represents the inner product, and i = 1,2, 3. It is noted that 


dh i 

5x7 ( x j- x k ) - 0 , 


(3.76) 


(3. 77) 


(3. 78) 


(3.79) 


(3.80) 
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where the i, j, k are cyclic indices 1, 2, 3. The slope of the projected limit cycle in the x k 
plane is 

dx k dh i /dx i 

dXj - dh^dx k ' (3.81) 

This is seen from the fact that 

dh. dtq dlq 

h, = TTT *i + 357 *j + dT k k k (3.82) 


in conjunction with Equation (3.80). Equation (3.81) corresponds to (3.66), and the second-order 
case reasoning is applicable. Of course it is true that in this two dimensional projection, the 
limit-cycle trajectory may cross itself (even though the limit-cycle trajectory in three- 
dimensional space does not) at points that are not system equilibrium points. That is, dh i /dx j 
= dh^dx^ = 0 implies an indeterminate planar slope at the point, but the system slope in three 
dimensions is parallel to the x £ axis. The essential idea is that the gradient vector in the x. x k 
plane is perpendicular to the projected trajectory, which in turn indicates the necessity for at 
least one peaking of the surface h i within the 
projected limit cycle. See, for example, Fig- 
ure 3.10. This Vh. = 0 peak corresponds to an 
tq = 0 equilibrium point as previously shown 
for the second-order case. Hence, at least 
one equilibrium point lies within the general 
limit cycle. 

Finally, the transition to the general 
n th -order system is made, where n > 2. The 
general form of the equation is 


projection c£ limit cycle 



*i “ f i < X I’ X 2’ * " ’ X n) 

i - 1, 2, ' • • , n . (3.83) 

Assume some initial conditions x £ (0) yield 
the limit cycle x. (t) = x £ (t +r) that is a 
closed bounded curve in n-dimensional state 
space. The projection cf the trajectory in 
the Xj Xj plane is represented by 

n— 1 rt 



n-2 


(3.84) 


Figure 3.10— hj(xj , x k ) = c s Curves. 



where the i’s are n cyclic variables 1, 2, * * • 5 n > and the c's are constant. There are n(n l)/2 
such expressions as in (3.84). Defining the field and gradient vectors as 


S 



and 


Vh. 

1 


1 1 2" ' 1 n- 2 



then 



(3.85) 


(3.86) 


(3.87) 


along the projected limit cycles. All cf the previous arguments hold true such that at least one 
^Fi = 0 point lies within all of the projected limit cycles. These points are equilibrium points 
the same as previously shown, since 


dx . 
fr 

dx. 


f . 


dh. . 

1 l 1 2” 1 n- 2/ 


d x . 


f . 


dh 


• • .. . • /dx. 

*1*2 1 n— 2/ *n 


Equation (3.88) is derived from (3.83) and the following two equations: 


(3.88) 


dh. 


<3h. 


dh : 


h. 


dx. 

1 1 


dx. 


dx. 


X; 


0 (3.89) 


and 


3h. 


' l n-2 


dx. 


K-,' *■„) = ° 


(3.90) 
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for j = ij , i 2 , * • • , i n _ 2 . Hence, the previous arguments from the second- and third-order 
cases still apply to the general n th -order case. 

It is therefore concluded that, for a general n th -order system that follows the ground rules 
of Chapter 2, if a limit cycle exists, then there is at least one equilibrium point that lies within 
the limit cycle. That is. 


x. 




< X 


< 


X 


iM 


i equib 


1 


1, 2, • • • n . 


( 3 . 91 ) 



4. EXAMPLES 


In this chapter, several examples are presented that illustrate the ideas of the previous 
chapters. Included also are "counter-examples" that apparently violate some of the theorems 
of Chapter 3. With each counter-example is an explanation of those ground rules in Chapter 2 
that are broken. 


Example 1 . First, a linear system is considered. The fourth-order differential equation 
is 


x ( t ) = -5<* 2 x(t) - 4w4 (x( t ) - x 3 ) , 

where w and x a are arbitrary real constants. Assume that initial conditions 

x(0) = x a + K 2 /4 

x(0) = 0 

'x(O) = 0 

'x(O) = 0 , 

where K is an arbitrary real constant. The solution is 

K 2 K 2 

x(t) X a + -j- cos wt - yj COS2CA , 

which is a limit cycle with period r = 2rr/oj. The periodic derivative functions are 



~K 2 03 

K 2 03 


i(t) 

= 3 sinwl 

+ 6 S1 

in 2cot 

x(t) 

-K 2 ^ 2 

K 2 CO 2 

-i- ^ 

cos 2oit 

— j coswt 

T 

x(t) 


2K 2 o> 3 


- “■ 3 ^ s i n w t 

_ 3 

sin 2oit 

x(t) 

K 2 w 4 

4K 2 « 4 

1 cos 2oJt 

= ”3 COSWt 



(4.1) 


(4.2) 


(4.3) 


(4.4) 
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Taking the averages of all these functions yields 


X 


avg 


X 


(i) 

avg 


/• t + 27 t/oj 
ij x(\)dX = x a 


Cti 

277 


/*t+27T/o) 

W J x(i> = o , 


where i 


1, 2, 3, 4. From Equation (4.1), the equilibrium point is at 


(4.5) 


(4.6) 

That is, the centroid or 
as. proved in Chapter 3. 
ure 4.1, illustrates how 
x-axis. 


average point (which lies within the limit cycle) is the equilibrium point 
Note that the projection cf the limit cycle on the x-axis, shown in Fig- 
the equilibrium point separates the positive and negative % regions of the 



Figure 4.1 -Linear-System Centroid. 


Example 2 . The second example illus- 
trates how one may easily be misled when 
solving for equilibrium points. Consider the 
periodic function 

x(t) = sint (4*7) 

2 his is the solution to the differential equation 

x 2 (t)+x 2 (t) = 1, (4-8) 


path the initial condition x(0) = 0. Now, x avg = 0, but it appears from Equation (4.8) that 
X equib = ± by set h n g *(*) = 0. If this is so, then these two equilibrium points lie on the limit 
cycle. See Figure 4.2. This contradicts the definition of an equilibrium point. Furthermore, 
x (t ) is also the solution to the linear second-order differential equation 


x(t) = -x(t) , (4-9) 

and it is shown in Chapter 2 that, for linear systems, x avg = x equib (as well as the fact that only 
one equilibrium point exists for a linear system). This is not the case when the equilibriums 
are calculated from Equation (4.8) which is the integral of the linear Equation (4.9). (Of course, 
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all is proper when Equation (4.9) is used to 
calculate the equilibrium points.) A ground 
rule cf Chapter 2 is violated by Equation 

(4.8) . Equilibrium points are defined as 

0 = f ( 0 , 0 , ■ • • , x e ) for the system 

~ f (x< n_1 >, X^H -2 ) , ••• , x, x) , (4.10) 

where n is a positive integer. But Equation 

(4.8) cannot be put into this form, such that 
every initial x yields a unique x. Either time 
is involved explicitly, or the system con- 
tains memory such that 



+ Vl - X 2 

for 

77 

0<t<~2 , 

(3 + 4k) 77- , 

2 - t 

(5 + 4k ) 77 
< 2 

Vl - X 2 

for 

(1 + 4k ) 77 
2 - 

. ( 3 + 4 k) 77 
t < 2 



(4.11) 


where k = 0, 1, 2, • • •. Misuse of relationships such as Equation (4.8) must be avoided. 


Example 3 . The well-studied nonlinear Van der Pol equation, 

x = e (l -x 2 ) x - x , e > 0 (4.12) 

has as a solution a limit cycle that encircles the origin in the x - x phase plane (Reference 11). 


See Figure 4.3(a). Substituting y - -x + e(x-x 



where x = y. Equation (4.13) likewise has a 
limit-cycle solution encircling the origin 
(see Figure 4.3b). In both cases the origin 
is an equilibrium point of the system. In the 
next chapter some further properties of 
Equations (4.12) and (4.13) are investigated. 

Example 4. This example illustrates 
two things: (1) equilibrium points may be 

difficult to find, and (2) an average point 
(centroid) is not necessarily an equilibrium 


/3) in (4.12) yields the equivalent system 



Figure 4.3— Van der Pol Equation Limit Cycle. 
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point for nonlinear systems. Consider the system 


-ofl (x + B) x < -X 

x = ^(k 2 x ) 1/3 |x| < X 


(j-o) 2 (x-B) x > X , 


( 4 . 14 ) 


where u and K are arbitrary positive real constants, and B = 20 kAj 3 >X> 0. Assume the initial 
conditions 


x( 0) _ B + A/ oc > 

x(0) = 0 , 


( 4 . 15 ) 


where A - 12/2K/C; 2 . This yields a limit cycle of period t = ( 3 n + 8 )/<•>, of which one cycle from 
05 t < r is described by 


x(t) = 


B + — cos ait 

Oj 


; H) J 


“B - ~ cos a) t 

Cl) 1 


-i) 


( ^ ?)' 


A 

B + ~ COS £U( t - t) 


0 < t < A 


A < t < j ~ A 


2 “ A < t <“2 + A 


+ A<t<T-A 


A < t < 


( 4 . 16 ) 


k(t) = 


-A sin wt 


-3K t - 


A sin ai (t 


3K [t - 


3 T 

4 I 


-A s i n a;( t - r ) 


0 < t < A 


A < t < ~2 - A 


r - r 

"2 - A < t < "2 + A 


+ A<t<T — A 


T - A < t < T 


( 4 . 17 ) 
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-A co cos cot 


■H) 

Aco COS U! - 

(•-*) 


A < t < ~2 ~ A 


2 A < t < 2 + A 


"2 + A < t < r - A 


“Act) COS Ctj(t-T) t - A £ t < r , 


where x(A> - x, and A - W 4<d. These functions are plotted in Figure 4.4 Note that x(t), and 
its first two derivatives are continuous for all t. This is forced to be true by constraining 


A 

B +— , 37 


3.33 % .3535 % .707 


A = 12 B = 20 j=(^j- + 2jY2 A = ^-i2 K = o = ^ 


/ 3.33 


x v ' \ A 

y 

® " . / 12.32 / \ 15.65 \ 

t/2 " A / -r/2 / \ r/2 + A V 


Figure 4.4 — Example4 Functions Versus Time. 



x, and x at t = 0, A, r/2-A, t /2 + A, t - A, and r by the equations 


B + ~ cos oA 

Ct/ 


AsinoA 


A&) eos «A - 



(4.19) 


respectively. To simplify the algebra, ^Ais arbitrarily chosen as "A - 377 / 4 , and A, B, and t are 
solved as functions of co and K. 

Equation (4. 14) apparently has three equilibrium points: 


^equib = 0, H . 


(4.20) 


Also, from Equation (4.16), one calculates the average of x(t) and finds 


avg 


l 

r 


f 


x( t ) dt 


0 • 


(4.21) 


(It is also true that x e - x e - x a - x a = 0, as expected.) Thus it appears that the origin is the 
centroid of the system, and is an equilibrium point. But this point lies on the limit cycle. The 
phase-plane trajectory is plotted in Figure 4.5. Its equations are, 


(x+B) 

2 a> 2 t x 2 

= A 2 

for 

X 

1 A 
1 

X 

x 3 = 

27 Kx 2 


for 

X 


i 3 - 

-27 Kx 2 


for 

X 

<0 J 

(x-B) 

2 O) 2 t X 2 

= A 2 

for 

X 

> X 


for 


x| < X 


(4.22) 


A defect in this example is that Equation (4.14 ) fails the uniqueness test of Chapter 2. That 
is. Equation (2.6 ) and (2.7) are satisfied, but the Lipschitz condition of Equation (2.8) is violated 
at the origin, since a Lipschitz constant must be greater than (x)~ 1/3 . In fact, x = 6(k 2 x) 1/3 
has three possible solutions from initial conditions at the origin. They are 

x 1 = 0, ±Kt 3 , (4.23) 
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of which only the first is an equilibrium solution. The last two are part of the limit cycle of 
Equation ( 4 . 16 ). 

The system is made unique if initial conditions that lie on the limit cycle shown in Figure 
4.5 (which includes the origin) are not considered. By varying A in the initial conditions of 
Equation (4.15), a family of limit cycles is obtained (see Figure 4 . 6 ) . The heavy figure-8 curve 
is the forbidden limit cycle where A = 12y2 YJ co 2 . However, it is not a separatrix (Reference 6 ) 
that separates regions and passes through equilibrium points so that it takes an infinite time to 
traverse. (A separatrix exists in the nonlinear pendulum example 8, and also example 11 , both 
of which are discussed later in this chapter.) 

If the initial conditions of Equation ( 4 . 15 ) are 


x(0) = + B ± A/co 

x( 0 ) = 0 , 


( 4 . 24 ) 


x 3 = 27K x 2 
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Figure 4.6— Phase -Plane Portrait for Example 4. 


such that A/co<B -x, then linear harmonic motion exists such that x(t) - +B+ (A/cu) cos The 
limit cycles encircle the ±b equilibrium points, which are also the centroids of the limit cycles. 
However, for the same initial conditions of (4.24)but B - X <A A,>< 12/2K / o? , the resulting limit 
cycles are still contained within the forbidden figure-8 limit cycle and still encircle the ±b 
equilibrium points; but these points are no longer their centroids. The centroids actually occur 
on the x-axis, but closer to the origin than ±B . For example, consider the X - 0 limiting limit 
cycle in the right-half plane. (The X - o curve is chosen even though it is part of the forbidden 
limit cycle, to avoid elliptic integrals.) x avg is calculated as 




(4.25) 


where 

12K^2 

A - ^ — 


B 


2QK 


377 

4 co 


, ^ 377 + 8 
T “ 2w ’ 


(4.26) 
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such that r' is half of the original forbidden 
limit cycle period r. Performing the inte- 
gration yields 


xavg 


B 


24 

1 5( 3tt + 8 ) 


0.725B , (4.27) 


which is not the equilibrium point, x e - B, 
inside the limit cycle. (When x avg and x a 
are calculated, they of course are zero.) 
This limiting limit cycle is shown in Figure 
4.7. 



Figure 4.7— Limiting Limit Cycle for Example 4. 


Example 5 . In Chapter 8, a mathematical on-off controller (of the sgn(x) type) is modified 
to match physical characteristics that have no jump discontinuities (see Figure 8.4). This ex- 
ample illustrates how a non-physical system may raise difficulties. Consider the system de- 
fined by 


x-K for x>0,0<x<K 

x - -k for x>0,-K<x<0 

x - -x elsewhere , 


(4.28) 


where K is a positive real constant. No equilibrium point exists, since x and x are never both 
zero simultaneously. However, many limit cycles exist, as shown in the phase-plane plot of 


Figure 4.8. This contradicts the theorem in 
Chapter 2, which states that if a limit cycle 
exists then an equilibrium point exists. The 
theorem is violated because the Fipschitz 
condition in Equation (2.8) is violated by the 
discontinuous jump of x. One cycle of x ver- 
sus time is plotted in Figure 4.9 for the 
heavily drawn limit cycle going through the 
origin in Figure 4.8. The jump discontinuity 
is re'medied by allowing x to vary anywhere 
between ±K (including zero) along the bound- 
aries of Equation (4.28), that is, x - 0 and 
x — 0. Thus, there are an infinite number 
of equilibrium points for x = 0 when 
i x | < K. (The origin lies inside all the limit 
cycles.) 



Figure 4.8— Phase Portrait for Example 5. 
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Example 6. This interesting example yields finite closed solution curves in the phase plane, 
but they are not limit cycles. The one point common to all the curves is the only equilibrium 
point of the system. The defining equations are: 


y - 2xy . 


The solution trajectories satisfy the relationship 


(4.29) 


x 2 + (y - c) 2 = c 2 ■ (4.30) 

Equation (4.30)represents a family of circles of radius c and with center on the Y axis at Y - c. 
From Equation (4.30) 


x 2 + y 2 
c 2y 

Substituting this expression for c in the derivative of (4.30), 

2xx + 2(y - c) y = q , 


yields 



2xy 


(4.31) 
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thus proving that Equation (^.SOjrepresents 
the integral solutions to ( 4 . 29 ) . The phase- 
plane curves cf ( 4 . 30 ) are shown in Figure 
4 . 10 . From Equation (4.29), the origin is the 
only equilibrium point. It takes an infinite 
time to traverse any of the closed circles. 
This example appears on page 10 of Refer- 
ence 8 , in Hahn's discussion of unstable 
equilibrium points. 

Example 7. Example 5 presents a sys- 
tem with a limit cycle but no equilibrium 
point. Example 6 presents a system with an 
equilibrium point but no limit cycle. Exam- 
ple 1 presents a system with both an equi- 
librium point and a limit cycle, and this ex- 
ample presents a system with neither. 
(Equation ( 3 . 18 ) in the previous chapter is 
also such an example.) The system under 
consideration is: 

x - xx - 2 . ( 4 . 32 ) 

Obviously, no finite real x exists such that 
"" = £ = 0. It is seen from Figure 4. 11 that 
any phase-plane trajectory crossesthe x - 0 
curve onlyonce, and the curve itself does not 
ever cross the x-axis. Hence a limit cycle 
does not exist. 

Example 8 . This example is the classic 
system of the simple pendulum, 

x" = -X. 2 sinx . ( 4 . 33 ) 

The second-order nonlinear system ( 4 . 33 ) 
yields an infinite number of isolated equi- 
librium points alongthe x -axis, i.e. x e = +imr , 
where n = 1 , 2 , • • •. Equation ( 4 . 33 ) can be 
integrated once, since 

.. . dx 

x = x dif 




Figure 4.1 1 -Zero-Slope Curve for Example 6. 


( 4 . 34 ) 
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Thus, 


x 2 - 2 k 2 cosx + K . (4.35) 

These solution trajectories are plotted in Figure 4.12. For initial conditions such that |k| <2A. 2 , 
limit cycles exist. K = -A 2 are the equilibrium points, and K = A 2 is the darkened separatrix 
curve which is the boundary between periodic and non-periodic solution curves. The limit cycles 
encircle a stable equilibrium point in between two unstable ones. Another interesting property 
of this nonlinear system is that the period of a limit cycle depends on the initial conditions (as 
opposed to the constant period in the linear Example 1). That is. 



Figure 4.12— Phase-Plane Portrait for Simple Pendulum. 


38 



Letting X (t 0 j - o , x(t 0 + r/4) - a, x (t 0 + r/4) - 0 (••• 0 - 2K 2 cos a + k) , and sin(x/2)-sin (a/2 ) sin 9, 
yields 


T 


= v4A 


f 72 dfl 

' 0 ^1 - sin 2 j sin 2 6 


( 4 . 37 ) 


which is a complete elliptic integral of the first kind, and is a function of a, which in turn is a 
function of K or the initial conditions. 


Example 9 . Example 8 contains an infinite number of isolated equilibrium points and limit 
cycles; Example 9 contains an infinite number of non-isolated equilibrium points, but no limit 
cycle. Consider the system, 


x = (x - y) (l - x 2 -y 2 ) 
y = (x + y) (l - x 2 - y 2 ) 

Choosing a Liapunov function (Reference 8, p. 18), 


( 4 . 38 ) 


. 2 = 


( 4 . 39 ) 


yields 


V = 2(1 - r 2 ) r 2 


( 4 . 40 ) 


Therefore, if 0< r 2 < 1 , then V > 0 and r 2 -l ; and if 1< r 2 , then v < 0 and r 2 - 1 as time increases. 
That is, the origin is an unstable equilibrium point, and anywhere on the unit circle in the x, y 
state plane are stable equilibrium points for the system ( 4 . 38 ) . The solutions to ( 4 . 38 ) in polar 
coordinates are 


r(t) - 


-1 1/2 


1 +Kj e~ 2t 


*(t) 


In 


K 2 

(l +Kj e” 2t ) 1/2 


( 4 . 41 ) 


where 


r - + Vx 2 + y 2 

9 = l an' ’I 


(4.42) 
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and 


K, 



K 


2 



(4.43) 


where r 0 and & 0 are the values of r and & at t - 0. Note that from Equation (4.41), 


lim r(t) 

t— W 


r Q ?0 


1 , 


and 


lim #(t) - In K, 

t-»C0 A 


There are no non-constant limit cycles in this system. Initial conditions on the unit circle give 
rise to constant solutions. A few trajectories are plotted in Figure 4.13. (Note that this system 
does not violate Chapter 2 ground rules, since every initial condition gives rise to a unique 
solution.) 

Example 10 . When an n th -order system is described by one n th -order equation such as 
Equation (2.9), then the equilibrium points are found by setting all the derivative terms equal 
to zero. When describing the system in the more general state variable form cf n first-order 



Figure 4.13— State-Space Portrait for Example 9. 


equations such as Equation (2.1), the n first- 
derivative terms are set equal to zero in 
order to find the equilibriums. It is not in 
general true that setting n arbitrary deriva- 
tive terms to zero (including second and 
higher derivatives cf any one cf the variables) 
for the state variable formulation, yields 
equilibriums. This example is an illustration 
of this fact for a second-order system. 

Consider the system described by 
. _ 

X I ■ X 2 4 T ( 1_X 1 2 ' X 2 2 ) 

(4.44) 

x 2 

X 2 = “ X 1 +_ 2 _ ( 1 _ X 1 2 _ X 2 2 ) ■ 
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This system has a limit cycle which is a circle of unit radius centered at the origin in the x t , 
x 2 state plane. That is, if 


Xl 2 (0) + x 2 2 (0) = 1, 


( 4 . 45 ) 


then 


x j ( t ) x 2 (0)sint + x j ( 0 ) cos t 
x 2 ( t ) = x 2 (0)cost- Xl (0)sint 


( 4 . 46 ) 


which is the unity-circle limit cycle with period r = 2 n, From Equation (4.44), there is only 
one equilibrium point, and it is at the origin (which of course, lies within the limit cycle). It 
is found by setting Xl - x 2 = 0, which yields only one solution, x t = x 2 = o. In fact, for any 
initial conditions other than at the origin, the system approaches the stable unity limit cycle 
with increasing time. That is. 


r(t) 



1 


+ 

f 1 > 

e-t 

(r 2 ( 0) ) 


( 4 . 47 ) 


0(t) = 0(0) - t , 

where r = [x 2 +x 2 2 ] 1/2 and 0 - tan -1 (x 2 /Xj), and 


1 i m r ( t ) 

t-*a> 


r (0 0 


1 


( 4 . 48 ) 


That the unity circle is a stable limit cycle is also shown by choosing a Liapunov function 


V ^ r 2 = X j 2 + X 2 2 

V - r 2 (l - r 2 ) , 


( 4 . 49 ) 


so that 


V > 0 for r 2 < l(r / 0) , .-. r 2 - 1 as t - co 


V < 0 for r 2 > 1 , 


/. 1 - r 2 as t - 20 


( 4 . 50 ) 


v = 0 for r 2 = 1 or 0, r remains 1, 0 for all t > o . 
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However, what happens if an “equilibrium” point is defined by x x - x x 
for the second-order system (4.44)? That is. 


0 or x 2 = x 2 = o 


and therefore, 


*i " f i ( x i- x 2 ) 

^2 f 2 ( x 1 * X 2) * 


a f i di 1 

X 1 <9Xj ^ 1 + (9 x 2 ^ 2 ~ Sx (xji X 2 ) 


3f 2 df 2 

x 2 _ dx x + <?x 2 ^ 2 - ®2 ( x l’ x 2 ) 


Hence, for x x = x x - 0 , 


f i ( x i> x 2 ) _ 0 

gl (xj, x 2 ) = 0 


(4.51) 


(4.52) 


(4.53) 


is satisfied. Thus, from Equation (4.52), if f 0 (since then f x - f , - 0, which is the equi- 
librium point already defined), then 


f i ( x i> x 2 ) - 0 

3fi 

dx 2 ( x i- x 2 ) ' 0 

For the system (4.44), the above equations (4.54 ) are 


Solving (4.55)yields 


0 = x 2 + 2 ( 1 _x i 2 ' x 2 2 ) 

0 = 1 - Xj x 2 


x i - iy 2 




+ Vs 


X 2 = l/x l 


(4.54) 


(4.55) 


(4.56) 
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However, the points given by Equation (4.56 ) yield non-zero x 2 and x 2 , that is. 


/5 

XI 


X - 


27*5 


(4.57) 


Hence, the points of (4.56)which satisfy x t = x x = 0, for the system (4.44 ) are not equilibrium 
points. (A similar situation occurs when x 2 = x 2 = 0 .) That is, setting x x = Xj = 0 in a sys- 
tem (4.51), does not necessarily yield an equilibrium point. 


In Figure 4.14, several curves are plotted in the x t , x 2 phase plane, f x (x x , x 2 ) = o and 
f 2 ( x i> x 2 ) ” 0 are the two S-shaped dashed curves. Note that they intersect only at the origin, 
which is the only true equilibrium state of the system. The S-shaped solid curve isx'j = g 1 (x 1 , x 2 )=0. 
It intersects the f 1 (x 1 ,x ? ) = o curve at two mirror-image points, x x = +1.27, x 2 = +0.786. To 
demonstrate that these two points are not 
equilibrium points, a system trajectory of 
(4.44 ) is shown with initial conditions 



(4.58) 

x 2 (0) = l/xj (0) = 0.786 

as the heavily drawn curve. The arrow- 
heads are in the direction of increasing time, 
and in several units of time the system has 
settled into its limit-cycle behavior. 

Example 11 . This last example illus- 
trates two things. One is the argument in 
Chapter 3 that the point where the gradient 
of j;he surface g(x) vanishes is an equilibrium 
point, and the other, like Example 4, is that 
the centroid of a limit cycle is not neces- 
sarily an equilibrium point for general non- 
linear systems. 

In the xy plane, consider the locus of a 
point P(x, y) such that the product of the dis- 
tances from P(x, y) to two points on the 



Figure 4.14— State-Space Curves for Example 10. 
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x-axis at x - ±a, is constant, say b 2 . In polar coordinates, the curve is expressed by (Refer- 
ence 12, Section 7.3) 


g(r , 9 ) ~ a 4 -b 4 + r 4 -2a 2 r 2 cos 10 - 0 , 


( 4 . 59 ) 


or, in x, Y coordinates, 

g(x, y) = a 4 b 4 + (x' + y 2 ) 2 - 2a 2 (x 2 -y 2 ) = 0 . (4.60) 

The graphs of the curves represented by Equations (4.59) and (4.60) for different values of the 
ratio b/a are shown in Figure 4.15. The curves in both (a)and (b) are called lemniscates, and 
the curve in (c) consists of two separate closed portions known as "ovals of Cassini." In Fig- 
ure 4.16, g(x, y) + b 4 is shown as a surface. The gradient vector as introduced in Chapter 3, is 

ry - if -v , § 3 . ~ 
v g - dx x + dy y 


- 4x(x 2 +y 2 - a 2 ) x + 4y(x 2 + y 2 + a 2 ) y . ( 4 . 61 ) 



Thus the gradient vanishes at three points, 


x - 0, +a 

y = 0 


( 4 . 62 ) 
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which correspondto the saddle point, and two 
minimum points of the surface shown in Fig- 
ure 4.16. 

The differential equations that uniquely 
gives rise to these curves are 


g(x, y) + b 4 


x(t) - y(t) (x 2 (t) +y 2 (t) + a 2 ) 
y(t) = -x(t) (x 2 (t) +y 2 (t) - a 2 ) , 


(4.63) 


since the solution curves are the limit- 
cycle trajectories 



Figure 4.16 — Surface of Example 11. 


(x 2 +y 2 ) 2 - 2a 2/ x 2 -y 2 ) = k 


.k> 0 leminiscate 
limit cycles 


k-0 

Separatrix 


-a 4 < X < 0 

Ovals of cassini 
limit cycles 


k=-a 4 

2 stable equilibrium points 



(The A = 0 curve is a separatrix for the different types cf limit cycles; it is not in itself a limit 
cycle since it passes through the equilibrium point at the origin.) The equilibrium points for 
Equations (4.63 ) are the same three points as given in Equations (4.62), that is, where Vg = 0 
there is an equilibrium point. The trajectories are plotted in Figure 4.17. Note that there is 
always an odd number, one or three, equilibrium points (i.e. surface peaks) within each limit 
cycle. 

Equations (4.63 ) were programmed on a digital computer, with a = 1, for a family of initial 
conditions. For each run. 


X avg 



x(t) dt 


(4.65) 


and t was computed, where r is the limit-cycle period. These values are listed in the Table in 
Figure 4.18. Note that both x avg and t are different for each limit cycle, and that 

X avg ^ X equib (4.66) 


x(0) 

T 

Xavg 

0 

00 

0 

±0.1 

6.72 

k0.468 

±0.2 

5.38 

±0.584 

±0,3 

4.66 

±0.675 

±0.4 

4.17 

k0.754 

±0.5 

3.82 

±0.822 

k0.6 

3.57 

±0.880 

±0.7 

3.38 

±0.929 

±0.8 

3.25 

k0.966 

±0.9 

3.17 

±0.991 

±1.0 

m 

±1.000 

±1.1 

3.18 

±0.989 

kl.2 

3.31 

±0.948 

±1.3 

3.67 

±0.856 

kl.4 

5.38 

±0.582 

± /2 

00 

— 


a = 1 
y(0) = 0 


avg 


= 0 
= 0 


Xavg 

y = 0 

J avg 

for all 
cases 


x(0) 

r 

Xavg 

±1.5 

6.38 

0 

±1.6 

4.58 

0 

kl.7 

3.60 

0 

±1.8 

2.96 

0 

±1.9 

2.50 

0 

±2.0 

2.16 

0 

±2.1 

1.89 

0 

±2.2 

1.67 

0 

±2.3 

1.49 

0 

±2.4 

1.34 

0 

±2.5 

1.21 

0 

±2.6 

1.10 

0 

k2.7 

1.01 

0 

±2.8 

0.92 

0 

±2.9 

0.85 

0 

±3.0 

0.79 

0 


Figure 4.1 8— Table of Computed Values for Example 11. 
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5. MISCELLANEOUS RELATED PROPERTIES 


This chapter is devoted to some miscellaneous discussions and concepts of ordinary dif- 
ferential equations and limit cycles that do not fit into the mainstream of the previous chapters. 
These discussions are included for their interesting properties, as an aid toward abetter un- 
derstanding of the overall topic, and possibly, to generate future ideas in this area. 


(1) For a second-order system, 

i(t) = f(i(t), x(t)) , (5.1) 

a solution trajectory in the x-versus-; phase plane that linearly (i.e., as a straight line) crosses 
the x-axis (i.e., x = 0) at a non-infinite slope, signifies that the crossing point is an equilibrium 
point. Such a crossing is illustrated in Figure 5.1. 



dx x _ f(x, x) , 
slope K - dx - ^ - x 

then for K finite at x - 0, f (0, x) is zero. Thus, this point is an equilibrium point. 

(2) For systems (obeying the ground rules of Chapter 2 ) of the form 

x^O - f , x( n-2) , • • • , x) , 


(5.3) 


(5.4) 
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only n 1.2 yields non-trivial limit cycles. That is when n - 0, Equation (5.4) is the algebraic 
equation x = f(x), which yields a constant solution with respect to time. 

ff n = 1, then 


x = f(x) . (5-5) 

Hence, when a solution x(t) is at a relative maximum or minimum value, x(t) is zero. But 
Equation (5.5) yields only unique solutions for any initial condition. Thus, x min = x Max = Xe i s 
an equilibrium point. That is. Equation (5.5) does not yield a time-varying solution with finite 
extremum values; hence, no finite limit cycle exists. 


(3) In general, it is very difficult to find limit- cycle solutions for nonlinear differential 
equations. This is true for even a second-order system, although the well-known Bendixson 
criterion (Reference 14) sometimes indicates where a limit cycle cannot exist. That is, if 


df i df 2 

7 ( x i- x 2 ) + 3)77 ( x i- x 2 ) 


dx 


(5.6) 


does not change its sign within a domain D of the x lf x 2 state plane, then no closed trajectories 
exist in that domain. A proof by contradiction is simple, since if a closed curve F exists in D, 
then by Gauss's theorem, for D' the domain bounded by F, 

dx 2 ~f 2 dxj) t 0 . (5.7) 


(5.8) 

Thus, fj dx 2 - f jdxj = 0 , which contradicts (5.7). 

For higher-order systems, Hahn refers to "rotating Liapunov functions" (Reference 8) of 
the form 


JI 


df t 3f 2 \ r 

3x7 + 5x7/ dx i dx 2 = J r ( f i 


But, along the path F, 


x i - f i ( x i> x 2 ) 
i 2 - f 2 (x,,x 2 ) 


u( x ) = f(x)/g(x) . (5.9) 

u(x) is interpreted as an angle-coordinate in a suitably chosen cylindrical coordinate 
systemunder certain assumptions about the pencil of surfaces determined by f(x) -cg(x) = 0, 
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i.e., by u - c - constant. As an example, for a third-order system, u - Xj/x 2 yields a pencil 
of surfaces which are planes in three-dimensional state space. If u> o, the corresponding mo- 
tion rotates about the axis of the pencil of surfaces always in the same rotational sense. Mo- 
tions always remaining within a torus-like domain infer the existence of almost periodic motion. 
The precise analytical formulation of rotating Liapunov functions is given in Reference 13. This 
concept is not necessary for two-dimensional phase-plane considerations, since the existence 
of periodic solutions is guaranteed if the phase trajectories remain in an annular domain. For 
second-order systems, the concept of "contact curves," as introduced by Poincare, is appealing. 

For the system (5.8), an arc A (open or closed at one or the other end) is said to be "with- 
out contact" wherever A. contains no equilibrium points and the (unique) solution trajectory path 
through any point of A. is never tangent to A at that point (Reference 11). Assuming these prop- 
erties hold for every point of a simple closed curve L, then L is called a "cycle without contact." 
This implies that once a solution trajectory enters (leaves) the region bounded by L, it cannot 
leave (enter) it. Furthermore, a structurally stable system has only a finite number of limit 
cycles in the region bounded by L (Reference 15). 

If a cycle (or arc) without contact is expressed by 

L^Xj, x 2 ) = constant , (5.10) 


then the tangent to L at any point is 


dh/dx , 

” dh/dx 2 ' 

For L to be without contact, the field vector 


(5.11) 


s = f ,*i + f 2* 2 ( 5 - 12 ) 

(see Chapter 3) is not tangent to L (and S is not equal to zero onL). However, contact curves 
C[x 1P x 2 )that are tangent to S are found by equating the slope of Equation (5.12) with that of 
c ( x i> x 2 ) given by (5.11) with L = C. That is, 

d C dC 

f 1 ( X l’ X 2 ) <9Xj ( X l* X 2 ) + ^2 ( X l’ X 2 ) dx 2 ( X l> * 2 ) ~ 0 ' (5.13) 

To illustrate this method of contact curves, a family of concentric circles centered at an 
equilibrium point, assumed to be at the origin, is considered. Thus, 

C = . x p + x 2 2 , (5.14) 
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and Equation ( 5 . 13 ) yields the locus of points, i.e., the contact curve, where the circles are 
tangent to the integral curves of ( 5 . 8 ). If a limit cycle exists at all, it lies in a ring domain, 
with center at the equilibrium point, whose boundaries are the innermost and outermost circles 
of radii r min and r MAX respectively, which touch the contact curve. (See Reference 6 , Figure 
2.2 for a numerical example.) Note that the equilibrium point, f 2 = f 2 = 0, is a solution to 
( 5 . 13 ). Thus, if it is considered part of the contact curve, r min = 0, and the limit cycle en- 
circles the equilibrium point as previously proved. The r MAX -radius contact circle curve is a 
measure of how large a region confines the limit cycle. If this contact-surface concept is ex- 
tended to higher than second-order systems, then limit cycles are contained in the hypervolumes 
within the largest norm hypercontact-surface. That is, this dissertation proves that limit cy- 
cles encircle equilibrium points so that only regions with equilibrium points need be considered. 
For practical systems, one expects the regions to be small. Contact surfaces mathematically 
define the boundaries cf these regions. 

Instead of the above topological approach, another viewpoint of circular contact curves is 
obtained dynamically by defining 


R(t) - |/xj a (t) + x 2 2 (t) ( 5 . 15 ) 

Then R = 0 when R(t) equals R min and R MAX (the minimum and maximum radii of the contact 
curve, respectively) and R min <R£R MAX . (A question that arises is: "Can any information about 
an 'average limit cycle', R avg , be obtained without knowing Xj (t)or x 2 (t)?") This approach 
sometimes aids in plotting limit cycles. 

For example, consider a second-order system of the form 


x = f(x, x) 


( 5 . 16 ) 


where R = fx 2 + x 2 . Then R = x(x + x)/R = 0 when 

x = 0 or x + x = 0 . ( 5 . 11 ) 

For the Van der Pol equation in the previous chapter (Equations ( 4 . 12 ) and ( 4 . 13 ) with 
y ~ x), x ~ 0 yields R = Q That is, the x-axis crossings of the limit cycle are inflection points 
of R(t). For Equation (4.12), with e - 1, x = (l -x 2 ) k-x, and 


R] 


X“ ± 1 


0 


Rl 


±i 


_ 2ex 3 

+ R 


( 5 . 18 ) 
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give the four alternating extremums of R(t) 
as shown in Figure 5.2(a). For Equation 
(4.13), with y - x, e - 1, x = x -x 3 /3 - x, and 

Kl x=±y3 = 0 

”i 2eV3x , c 

= ± -^~ ( 5 - 19 ) 

give the four alternating extremums of R(t) 
as shown in Figure 5.2(b). Also plotted in 
Figure 5.2 are the circles of radius R avg . 

(4) A great deal is known about the be- 
havior of phase portraits around singular 
points for second-order systems. (See Ref- 
erence 16, for example.) For higher-order 
systems very little is known. Among others, 
Poincare" (Reference 17) studied third-order 
singularities, Mendelson (Reference 18) in- 
vestigated special higher-order systems, and 
Haas (Reference 19) discusses complex sin- 
gular points. However, for most practical, 
complicated, high-order systems a computer 
is needed; also, engineering ingenuity plays 
an important role. (For example, time is run 
backwards and forwards to determine stabil- 
ity of equilibrium points, limit cycles, and 
alternating stable- unstable limit cycles.) 
This dissertation is intended to bridge part 
of the gap that exists between mathematical 
theory and engineering applications. Much 
more knowledge is needed. For example, 




(b) 


Figure 5.2— Properties of Van der Pd Limit Cycle. 


what are the necessary conditions for equilibrium points and limit cycles to exist? (Chapter 3 


showed that the system equations must allow x. to take on values of both polarities, since x. = o 

avg 

for periodic x.) Another interesting question is: When is a centroid an equilibrium point? 

(That is, when does f(x avg ) = 0, as proven for linear systems?) What about non-autonomous 


systems, and systems with time-varying parameters? How can equilibrium points be moved 
around in state space, that is, what are the sensitivity functions dn i Jq (parameter)? These, 
and many more questions remain to be answered. 
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Part II 


APPLICATION TO ATTITUDE CONTROL SYSTEMS 

6. INTRODUCTION TO ATTITUDE CONTROL SYSTEMS 

Attention is now focused on spacecraft attitude control systems. Of special interest are 
three-axis stabilized, either vertically earth-oriented or inertially oriented, earth-orbiting 
satellites. Normally, six positional degrees of freedom are associated with rigid-body motion; 
three translational and three rotational. The guidance and control requirements for navigating 
the center of gravity of a spacecraft with respect to a desired trajectory or orbit (or midcourse 
maneuvering or orbital station keeping) is not discussed here, although the techniques of this 
dissertation are definitely applicable. Instead, attention is restricted to the control of the ro- 
tational motion of a rigid body, or system of coupled rigid bodies, by means of devices that 
store or change the momentum of the bodies. Only the three rotational-position degrees of 
freedom per rigid body need then be considered, and since the rigid bodies may have sources 
of torque acting upon them, there are three rotational-rate degrees cf freedom yielding a six- 
dimensional system per rigid body, with constraints relating the bodies. (If a torque applied to 
a single rigid body is not a pure couple, then the trajectory of the center of mass changes and, 
in general, one must consider a 12-dimensional system. However, the torques and forces as- 
sociated with attitude control are usually extremely small compared to any that significantly 
change the trajectory cf a spacecraft. Asa matter of fact, one is often forced to cope with the 
small components of huge orbit stationkeeping thrust vectors (such as to keep a synchronous 
24-hour orbital period earth satellite in a proper orbit) which do not pass through the space- 
craft center of gravity, as attitude disturbance torques). 

It is emphasized at this point that the intention is to present not an all-encompassing de- 
velopment of attitude control and related background areas, but rather a broad framework that 
leads to system equations to which the equilibrium technique is applicable. The development is 
kept fairly general and practical, in order to understand the equations in a physical sense. A 
complete general attitude control system is shown in Figure 6.1. The system is divided into a 
designable spacecraft and its governing laws of motion. The first step taken by the attitude 
control engineer is to develop the system equations in a mathematically convenient form. This 
first step is discussed in the next chapter. In this chapter, the philosophy of a design-analysis- 
redesign-reanalysis-etc. approach is presented, for the sake of a physical understanding of the 
equations. Ideally, a paper design is desired that results in an optimum or near-optimum con- 
trol system that satisfies all the various spacecraft requirements. Once the paper design is 
completed, a synthesis-type problem would exist. Unfortunately however, owing to practical 
limitations and the lack of analytical tools for highly nonlinear systems, practical designs are 
based on engineering intuition and experience. The effectiveness of these designs demands the 
solution cf an analysis-type problem for which a computer serves as the main nonlinear tool. 
This important role of the computer is pursued further after a brief discussion of the nonlinear 
dynamics and kinematics of the system shown in Figure 6.1. 
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Figure 6.1 —Attitude Control System 

The dynamic behavior of a system is determined by the equations which relate the external 
torques and the angular momentum of the system. The momentum consists cf two components: 
the momentum cf the rigid body due to its angular rate W with respect to inertial space, and the 
stored momentum from internal rotating members. The dynamic behavior cf attitude-controlled 
systems is determined by three scalar nonlinear, first-order, time-invariant, ordinary differ- 
ential equations. These equations are derived in Chapter 7, and their scalar form, Equation ( 7 . 45 ) 
is 

I,Wi = (Ij-I 3 )W 4 W 3 +N, -h, +h 2 w 3 -h 3 w a 

lA - (l 3 -Ii)W 3 Wi rN.-h.+hjW, -h ( W 3 (61 ) 

I 3 w 3 - (irl^WjW, + N 3 - h 3 + h,W 2 - h 2 w, , 

where the symbols and assumptions are given in Chapter 7. These equations are sometimes re- 
ferred to as the Euler equations of motion and are heavily crosscoupled with nonlinear terms. 

As opposed to the well-defined dynamics equations, the formulation cf the kinematics of a 
system is more arbitrary in that there are many sets cf position variables to choose from in 
order to describe the attitude of a rigid body. The kinematics are relationships between the 
rates and one of these sets. In translational motion, the kinematics are simply 

— v . i — T 2, 3 r 

where the set is a 3 -parameter rectilinear coordinate system, For rotational motion, the kin- 
ematics become much more complex nonlinear coupled relationships, but still, only three in- 
dependent coordinates are necessary to specify the attitude. The three most commonly used 
sets of coordinates are now discussed. 
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The first set is the three Euler angles taken in any advantageous (but independent) sequence, 
such as azimuth, elevation, and "spin" cf a gun mount. See Reference 20' for example. The main 
advantage cf this set is that there are only three first-order differential equations relating the 
angles and the angular rates. They are very easy to visualize for reasonably limited angular 
motion (such as the roll, pitch, and yaw of a ship), and for sufficiently small angles become 
linearly decoupled, thus facilitating small-angle fine-pointing control analysis. The main dis- 
advantages are the non-uniqueness of the formulation, and the existence of singular points for 
large angles (90 degrees). Since computers cannot maintain accuracy for infinite quantities, 
Euler angles are not in general suitable for all-attitude simulations. As an example, taking a 
yaw, roll, pitch, (i.e., ^ , 4 , 6 ) sequence yields: 

W - (ty cos 6 - W, sin#)/cos<£ 

& - W 2 + Wj tan<£sin0 - W 3 tan<£cos 9 2 ) 

4> ~ WjCos5 + W, sine 

The second set cf coordinates is the nine direction cosines, (i, j = 1 , 2 , 3 ). The formu- 
lation contains nine defining first order differential equations which relate the positional direc- 
tion cosines to angular rates, plus six algebraic equations of constraint, thus yielding the three 
necessary independent position variables. Hence, there, are many equations and many variables, 
but they are all well-defined and well-behaved for all attitudes. It is easy to mentally visualize 
the physical meaning cf the direction cosine which enhances intuitive-type manipulation cf the 
equations, and lends itself to the functional developments used later in this dissertation. See for 
example Reference 21 . However, the direction cosine formulation requires unnecessary com- 
puter equipment for direct simulation. The nine differential equations take the matrix form, 




a n 

a 1 2 

a l^ 


0 

W 3 

-w 2 

d r - 

a 








dt [A] 

- dt 

a 2 1 

a 2 2 

a 23 


~ W 3 

0 

"1 



a 31 

a 32 

a 33 


_w 2 

-w, 

0 


where the a... are the nine coordinates. The 6 constraint equations are 


[A] [A] T 


10 0 
0 10 
0 0 1 


(6.3) 


(6.4) 


where [A] T is the transpose of [A] . 

The third set of coordinates is the four Euler symmetric parameters (Reference 22 ) whose 
formulation is a good compromise between the efficient but ill-behaved Euler angles, and the 
non-efficient but easy to use direction cosines. With only four defining first-order differential 
equations and one algebraic equation of constraint, the four variables and their derivatives are 



conveniently bounded and uniquely defined for all attitudes. They also lend themselves to a 
simple physical interpretation. This formulation is usually the most efficient to use in all at- 
titude analog simulations and then, if desired, the direction cosines are formed algebraically 
from the four symmetric parameters (Reference 1 23) . This is usually done to facilitate the 
sensors' simulation since the mathematical models of attitude sensors are most simply de- 
scribed via direction cosines. The four differential equations take the form 



where the equation of constraint is 



= 1 , 


(6.5) 


( 6 . 6 ) 


and A„ k 2 , \ 3 ,X. 4 , are the four position coordinates. (See Reference 23 for relationships be- 
tween all three sets of position coordinates.) 

There are additional sets of position coordinates which are used less frequently for de- 
scribing the kinematics, such as the Gibbs vector, Quarternions, and the Caley-Klein param- 
eters (Reference 24). They are sometimes used to advantage in specialized formulation of the 
control laws for paper-and-pencil analysis techniques. 

Thus the dynamics and kinematics are determined by well-defined mathematical relation- 
ships which close the loop around any rigid body as shown in Figure 6.1. It is the spacecraft 
hardware which is designed by the attitude control engineer to make the overall system work 
properly. The sensors (such as gyros, star trackers, sun sensors, earth horizon scanners, 
etc.) and the muscles (such as torque-generating pneumatic subsystems, momentum- storing 
flywheels, etc.) are usually selected very early in the design since they are the long-lead 
items. In fact, it is usually not until after the sensing and torquing components are being de- 
veloped that the overall simulation becomes more realistic, since then the hardware charac- 
teristics are measured. Hence, upon analyzing the system, the only functional design change 
which is comparatively easy and not too time-consuming is in the control laws. Modification 
of the sensors and muscles is possible to an extent, but the really flexible part of the system 
is the control logic. Therefore, after a system goes through its basic design, its performance 
is analyzed and redesign initiated in the event that specifications are not met. This cycle is 
continued until a satisfactory design is achieved. The computer simulation is essential to this 
type of design-analysis cycle, since, once the dynamics, kinematics, and sensors are simulated 
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(the most difficult and component-consumingaspects of the simulation), it is not difficult to 
change the simulation of the control laws. It is demonstrated in this dissertation how the com- 
puter can be intelligently and fruitfully used to find the potential trouble areas via system- 
equilibrium points, and further indicate how they can be remedied. 
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7. DEVELOPMENT OF ATTITUDE CONTROL EQUATIONS 


In this chapter, the necessary dynamic and kinematic equations are developed using the 
direction cosines which relate a body-fixed spacecraft reference coordinate system to an 
inertially-fixed reference coordinate system. The more general case of an "in-between" non- 
inertial orbital reference frame (for example an earth-oriented orbital frame which rotates 
one revolution per orbit in inertial space) is considered first. That is, for perfect attitude, 
the control system aligns the spacecraft reference frame to the orbital reference frame. The 
assumption is made that the plane of the satellite's orbit is inertially fixed. This is usually a 
good approximation since, during the time interval in which an attitude control system responds, 
the orbit precesses a negligible amount. Also, the 0.04 degree/hr earth rate around the sun 
and similar motions are ignored. If it is necessary to take long-term disturbance effects into 
account, a fourth reference frame is fixed (for example to the sun) and considered to move in 
inertial space (for instance with respect to a star). 

Since this dissertation is not intended as an all-encompassing text on attitude control and 
related areas, several additional assumptions are made, any of which can be readily removed 
by generalizing-type modifications to the equations. Orbits are taken as circular about per- 
fectly homogeneous spherical earths. No translational, but only rotational, motions are con- 
sidered with all reference-frame origins continuously coincident. The spacecraft is a single 
rigid body, and all inertial momentum storage devices have negligible inertias. Cross product 
of inertia terms, and sensor-actuator misalignment terms are omitted. 

For present purposes, the following three orthogonal right-hand Euclidean reference 
frames are defined. (See Figure 7.1). As previously stated, the only interest is rotational be- 
havior; thus all frames are assumed to have origins at the spacecraft's center of mass. The 
symbol " represents a spacial unit vector. 

The first reference frame is the inertial frame, z t , z 2 , z 3 , such that the z 1} z 3 plane de- 
fines the spacecraft orbital plane which was assumed inertially fixed. At some point in the 
orbit the positive z 1 axis points in the same direction as the velocity vector of the spacecraft's 
center of gravity, and the positive z 3 axis points to the center of the earth (local vertical). A 
quarter of an orbital period later, £ t points away from the center of the earth and z 3 is along 
the spacecraft velocity vector since the inertial frame does not rotate as the spacecraft goes 
around the earth. The z 2 axis is always perpendicular to the orbital plane throughout the 
orbital period and is to the right as one looks down the positive z 2 axis, i.e., z 2 = z 3 x z 1 . 

Hence the momentum cf the spacecraft due to its center of mass being in orbit around the 
earth is along the negative z 2 axis. 

The second reference frame is the orbital reference frame, y 1? y 2 , y 3 , such that: the 
positive y 3 axis always points along the spacecraft velocity vector in the direction of forward 
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motion; the positive y 3 axis points along the local vertical towards the center of the earth; and 
■the y 2 axis is coincident with the inertial z 2 axis. Thus, if aft) = n Q t, where ft 0 is the constant 
orbital rate in rad/sec, then a is the orbital angle defining the position of the spacecraft in its 
orbit. When the orbital reference and inertial frames are coincident, a = 0. The inertial and 
orbital reference frames rotate with respect to each other about the colinear y 2 and z 2 axes at 
a rate fi 0 . They are related by the matrix equation 


M 


cos a 0 sin a 


f 2 \ 

y 2 


0 10 


^2 

Wi 


- sin a 0 cos a 




or, in shorthand notation, 

(y) = [a] (z) . (7.2) 

The third reference frame is the spacecraft frame x y x 2 , x 3 . It is coincident with the 
principal axes of the spacecraft. When the spacecraft is perfectly controlled, the principal axes 
are aligned with the orbital reference frame. x 2 , x 2 , x 3 are the roll, pitch, and yaw axes, re- 
spectively. Because they are along the principal axes of the spacecraft (which is considered a 
rigid body), the cross product of inertia terms vanish. Therefore, the inertial matrix associated 
with the inertia tensor (Reference 21) becomes diagonalized with only three non-vanishing terms: 
I 1} I„ I 3 . These are the roll, pitch, and yaw principal moments cf inertia and are positive 
quantities. At any time, each of the body-fixed unit vectors may be resolved into components 
along the orbital reference directions. This fundamental concept allows one to express mathe- 
matically each of the spacecraft fixed unit vectors as a linear combination of the orbitally fixed 
unit vectors as 


3 



i = 1 


(7.3) 


where the coefficients a i . are time-varying if the spacecraft rotates with respect to the orbital 
frame. In matrix notation 


j*A 


3 11 ai 2 a 1 3 



5 2 

- 

a 21 a 22 a 2 3 


A 

y 2 

{*>} 


a 31 a 32 a 33 


[ysj 


or, in shorthand notation, 


(x) - [A] (y) . 

Since the x, y, and z coordinate systems are orthogonal, 


X. • X, 


6 .. 

ij 


y . • v • — 5 . . 

z . • z . = 5 . . 

1 3 


i, j - 1, 2, 3, 


(7.4) 


(7.5) 


(7.6) 
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where the Kronecker delta is defined as 


fl for i - j 
§ii to for i f j 


(7.7) 


Hence, multiplying both sides of Equation 73 by y k , yields 


3 ^ 

e • y k = J2 a o (?> ■' y «) = Z] * i > 


3 ., 

jk 


ik 


(7.8) 


j = 1 i = i 

Since the geometric definition of the inner product of two vectors is 


y k 


y k 


cos ^ 


= cos,; 


(7.9) 


vk 

where | | represents the magnitude cf the vector, and where is the plane angle between the 
two vectors, then 


a, = cos); (i> k = 1, 2, 3) (7.10) 

The a ik are the so-called direction cosines. This geometric interpretation is shown graphically 
in Figure 7.1 fork = 3. 

The transformation matrix [A] has many interesting properties from the point cf view cf man- 
ipulation and computation. Forming the inner product (xj ■ x k ) from Equation (7.3), one obtains 


3 3 

«k = 21 ^ 12 a,tm ^ 

j = 1 m= 1 


3 3 


EE- 

j = 1 m- 1 

Ee- 

j- 1 m= 1 


* km y i ‘ y » 


. . a, 3 . 

* * km j m 


■ E - ■*. 

j=i 

But from Equation (7.6), Xj • x k = S ik . Hence, 

3 

7~! a ii a kj = s ik . 

where i, k = 1, 2, 3. 


j = i 


(7.11) 


(7.12) 
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Equation (7. 12) represents six properties of [A] which will later be referred to as six 
equations cf constant. Equation (7.12) establishes that the sum of the squares of the elements 
of each row of [A] is unity, and the sum cf products of the respective elements of each pair of 



rows of [A] vanishes identically. When the six equations, (7.12), are written out, they are 



-.1 

+ aj 2 2 + 

3 1 2 3 

11 


a 2 2 l 

+ a 2^ + 

3 23 

= 1 


2 

a 3 1 

+ a 2 + 
3 2 

a 3*3 

= 1 


a 21 + 

3 1 2 a 22 

+ a 1 3 

a 23 

a n 

a 31 + 

a 1 2 3 3 2 

+ a 13 

a 33 

a 21 

a 31 + 

a 2 2 a 32 

+ S 23 

a 33 


The transpose of the matrix [A] is 



a n 

3 2 1 

a 31 

( 1 

> 

it 

a 1 2 

a 2 2 

a 32 


a 1 3 

a 23 

3 3 3 


(7.13a) 


(7,13b) 


(7.14) 


Clearly, the element in the i th row and j th column cf the product matrix [A] [AI T is 


Therefore, 




a u 

a 1 2 

a 1 3 


a ll 

a 21 

a 3 1 

[A] [AD T - 

a 2 1 

a 2 2 

a 23 


a 1 2 

a 22 

a 32 


a 3l 

a 32 

a 33 


a i3 

a 23 

3 3 3 


1 

0 

0 


0 

1 

0 


0 

0 

1 


= lul , 


(7.15) 


where [u] is the three-by-three identity matrix. Thus, [A] is a unitary or orthogonal matrix. That 

1S ' (AT 1 = [A] T (7.16) 

where [A] " 1 is the inverse matrix of [A] . Also, 

DetlA] = | A[ - +1 (7.17) 

since all reference frames are right-handed sets. 
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Another useful property of [A] is that 


a ij 1 cof ( a ij) > ( 7 - 18 ) 

where cof( a^) denotes the cofactor cf element a.. and is defined as the product of (-l) i + i and 
the determinant cf the square array formed by removing the row and column in which a., ap- 
pears from [A] . 2his is readily shown from 


[A ]' 1 


Adj [A] 
Det [A] ’ 


(7.19) 


where Adj [A] is the adjoint matrix of [A] formed by replacing each element a., of [A] by its 
cofactor and transposing. Substituting Equations (7.16) and (7.17) in Equation (7.19 ) yields 


[A] T = Adj [A] , (7.20) 

which is the same as Equation (7.18) . Writing out these nine relationships leads to 


a ll 


a 22 

a 33 

a 2 3 

a 3 2 

a l2 

.= 

a 2 3 

3 3 1 

a 2 1 

a 3 3 

3 1 3 

= 

a 2 1 

3 3 2 

a 22 

a 3X 

3 2 1 

= 

3 1 3 

a 3 2 

- a j 2 

a 3 3 

®22 

= 

a ll 

a 3 3 

■ 3 1 3 

a 31 

3 2 3 

= 

a 1 2 

a 31 

“ a ll 

3 3 2 

a 3 1 


a 1 2 

3 23 

" a 13 

3 2 2 

3 3 2 

= 

3 1 3 

a 2 1 

~ a ll 

3 2 3 

a 33 

■* = 

3 1 1 

a 22 

a 1 2 

a 2 1 


(7.21) 


The last useful property noted is another form of the six constraint Equations, (7.13). These 
follow directly from equating the elements of 

[A] T [Al = [U] , (7.22) 


which yields 



(7.23) 
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or, 


a ft 

+ a 2~l + a 3 2 l ' 

1 


3 1 2 2 

+ a 22 + a 32 

1 


a Vi 

2 , 2 — 

+ a 23 + a 33 

1 


a l 2 

+ a 21 a 22 + a 31 

a 32 

= 0 

a l3 

+ a 21 a 23 + a 3 1 

®33 

= 0 


a i2 a 13 + a 2 2 a 2 3 + a 32 a 33 = 0 . 


(7.24a) 


(7.24b) 


Many more interesting relationships can be derived (Reference 24), but for present purposes, 
the above suffices. 

It is often useful to have available differential equations, rather than algebraic equations 
of constraint, for manipulative purposes. From Equations (7.12) and (7.15), the following are 
obtained: 



i, k = 1, 2, 3 


(7.2V) 


or its equivalent 


[A] ft [A] t + |[« • [Al T = 0 . 


1 From Equation (7.17), | A| is constant. Therefore 



(7.15') 


(7.17') 


and [A] is singular, possessing no inverse.* From Equation (7.18), 


d 

a ii - dt [ cof a ij] . i J = 1,2,3. (7.18') 


*In general, Detd/dt[A] ^ d/dt Det[A], but for [a] orthogonal, [A] [a] T - [u], [A] [A]"*" + [a] [a1^" = 0, ■■ [A] 61 [a]"*" [a] . 

Thus, |A| = (-if |a| |A T | |a|= - |d/dt A T I, since the determinant of the product of two square matrices is the product of the de- 
terminants (Reference 25), n = 3, and |a| = 1 — But, the determinant of a matrix is equal to the determinant of its transpose, hence 
|A| - |A T | = |d/dt A^ | = - |A|. Therefore, |A| = 0. 
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Finally, from Equation (7.23), 


3 

(a. . a ik + a. . a lk ) =0 k, j = 1, 2, 3 , 

1=1 

which, when written out, yields 

'll a ll + a 21 a 21 + '31 '31 ® 

'12 a 1 2 + a 22 a 22 + a 32 '32 ® 

a 1 3 3 1 3 + ’ 23'23 + ' 33'3 3 ® 


a ll 

a 1 2 

+ 

a 2 1 a 2 2 

+ 

’31 a 32 

+ 

a u 

•12 + 

a 2 1 

-22 

+ 

a 31 

^32 

0 

a u 

3 1 3 

+ 

a 2 1 3 2 3 

+ 

a 3 1 a 33 

-1 

a ll 

4 13 + 

3 2 1 
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a 3 1 

^33 

= 0 

’12 

3 1 3 

+ 

a 2 2 a 2 3 

+ 

'32 a 33 

+ 

a 1 2 

’13 + 

3 2 2 

’23 

+ 

a 32 

’33 

0 


(7.23 1 ) 


(7.2 5a) 


(7.25b) 


The main point is, as „uted in Chapter 6, that there are only three independent generalized 
coordinates required for relating the rotation of one reference frame to another. By using 
direction cosines, there are nine elements introduced. Hence, out of all the foregoing relation- 
ships from Equation (7.12) on, there are only 6 independent constraining scalar relationships, 
but any of the algebraic or differential ones may be used for conveniently simplifying expressions 

The dynamics and kinematics are now developed, and the general equations of motion derived 
To begin with, 


m = cq Xj + o> 2 x 2 + CtJ 3 x 3 (7.26) 

is defined as the vector representing the rate of rotation of the spacecraft about its principal 
axes of inertia with respect to the orbital reference frame y . Since the orbital frame y rotates 
inertially about the y 2 axis at a rate -fi 0 , 

W = w - n 0 y 2 (7.27) 

where W is the inertial rotation of the body frame. From Equation (7.5), 

y = [Al T x . (7.28) 


Thus, 


+ a 22 X 2 


+ a. 
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hence, 


W 


{“l ^o a 12)^l + ( W ’ ^0 a 22) X 2 + (^3 a 3 2) X 3 

W x 1 , + W 2 + W 3 x 3 . 


('7.25; 


Qf course, if the orbital reference frame to which the attitude control system aligns the space- 
craft is inertial rather than earth-oriented, then n 0 = 0 and co = W. Or, for attitude motion 
during a time interval that is small with respect to an orbital period 27t/& 0 , so that K| >>fi 0 5 
orbital rate is neglected and a> % W. (Remember, |a.. j <1, i, j = 1, 2, 3.) 

>|f 

The total angular momentum, H, of the spacecraft consists of two parts. One is H s/c that 
is due to the rigid body rotating in inertial space, i.e., 



T $ <? T 

Xj x x 

0 

0 


”wi x 7 

H s/c = 4* • W = 

0 

I 2 £ 2 S 2 t 

0 


W 2 i 2 


0 

0 

T A /v T 
I 3 X 3 X 3_ 


5 \ 


= Il W l X l + IjMj + I 3 W 3 X 3 - ( 7 - 30 ) 

where 0 is the previously mentioned inertia tensor in dyadic matrix form. The second is h due 
to the momentum cf internal rotating parts such as reaction flywheels, rotating solar array 
paddles, tape recorders, etc. That is, 

H = H s/e t h ^ H, Xj t H 2 x 2 + H 3 x 3 . (7.31) 

Calling N the total external torque on the spacecraft, one can write 


N - Nj Xj + N 2 x 2 + N 3 x 3 . 


(7.32) 


*s/c is an abbreviation for spacecraft. 
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H and N can also be resolved into inertial coordinates via [A] and [a] . That is, if H and N are 
expressed as column matrices in the body frame, 


(H) = 


(N) = 


/HA 

h 2 

w 

K\ 

n 2 

\ N 3/ 


and as column matrices in the inertial frame, 


then they are related by 


(H ) 1 = 


(N) 1 = 



(H) = 

[A] [a] (H) 1 = 

[B] (H) 1 

(n) = 

[A] [a] (N) 1 = 

[B] (N) 1 


(7.33a) 


(7.33b) 


(7.34) 


where [Bl = [A] [a] . The elements of [B] , b.. , all have the same properties as , since [B] , 
like [A] , is an orthogonal unitary transformation mattix. 

Just as Newton's inertial law equating force and time-rate change of translational momentum 
is basic to linear motion, so is the fundamental equation of rigid-body motion. It is simply writ- 
ten in matrix notation as 


3t (H) 1 = (N) 1 


( 7 . 35 ) 


and is basic to all rotational motion. 
(Reference 21) as 


Equation (7.35) is written with respect to body coordinates 


d - 

JH + WxH = N , (7.36) 
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where 


(w) = 


K\ 

W 2 


\ W 3 / 


Equation (7. 36) is equivalent t o 


5E (h) + [p w ] Ch) = (n) , 


(7.37) 


(7.38) 


where 


[ p w] 


o -w 3 w 2 

W 3 0 -Wj 

-W 2 Wj 0 


(7.39) 


Note that |p w | = 0. Hence, [p w ] is a singular matrix and, as such, possesses no inverse. 
Also note that 


t P w] T = ‘[Pw] • 


(7.40) 


Now from Equation (7.35 ) and (7.34), 

(H) 1 = (N) 1 = [B] “ 1 (n) . 


(7.41) 


Also 


5t (h) 1 = {[Bl - 1 (H)} 

= £ {M' 1 } ® + [B]- 1 £ (H) • (7.42) 

Equating the right-hand sides of Equations (7.41 ) and (7.42), premultiplying both by [Bl , and 
realizing [B] -1 = [b] t , yields: 


gp (H) + [Bl ^ { [B] t } (h) = (N) . 


(7.43) 


67 



Equations (7.36), (7.38 ) and (7. 43) are all equivalent ways of stating the well-known three 
scalar Euler equations where 


(H) 


Mi + h i\ 
I 2 W 2 +h 2 
\ I 3W 3 + h 3j 


Written out, their familiar form is 

Ni = + (I 3 -I a ) W 2 W, + h, + h 3 W 2 - h 2 W 3 

N 2 = I 2 w 2 + (Ii-I,)W 3 Wi + h 2 + h, W, - h 3 W, 

t N f = ^ + (I 2 - Ii) W i w 2 . +>3 + h 2 W 1 - h t W 2 

V v V 

Rigid body Internal rotating parts 


It is clear from Equations (7.38) and (7.43)that 


CB3 dT [BJ T = [p w ] 


Premultiplying both sides by [B] ~ 1 yields 


dt [Bl T = [B]“‘[Pj 


Taking the transpose of both sides and using the identities 


[B] 


1 - 


[B] ' 


t P w] T = “[Pw] 


£ {ra T ) 

{[x] [y:} t 


IF [B] 


Cy] t [x] t 


yields 

^[B] = -[P W ][BJ.* 

‘Thus |d/dt [B]|= " |P W | |B I = - |P W |, and both [b] and [P w ] are singular matrices, as previously stated. 


(7.44) 


(7.45) 


(7.46) 


(7.47) 



The equivalent vector form of (7.47) is 


^ = b. x W , 


i - 1, 2, 3 


where 


b. 


Ti x i 


b 2i *2 


+ b 3i *3 


(7.48) 


( 7 . 49 ) 


Equation (7.47) or (7.48) represents nine scalar first-order differential equations. These nine, 
along with the three Euler equations, plus all necessary constraints, comprise the dynamics 
and kinematics of the spacecraft. They must now be put in the form used in Chapter 8. Firstly, 
it is assumed that for all practical purposes, any significant motion takes place in a sufficiently 
short time with large 5 , so that [A] = [B) and is ignored. Hence W = T, Also, it is assumed 
that any external torque generator, such as a pneumatic or mass-explusion system, is much 
larger in torque than the torques produced by the internal rotating elements. Hence, h = h = 0. 
One must be careful with this assumption since even for small h, ifxh can be large if W is large. 
A detailed computer simulation is usually necessary to show that the h terms do indeed cause 
only small perturbations to the system motion. This is most often the case for active systems, 
but caution must be exercised for passive or semi-passive systems. An additional assumption 
is to consider 2-axis positional control that aligns a principal axis to an orbital (inertial) axis, 
and controls the spacecraft rate about this axis, but not the angle. Two-axis control of this 
nature is all that is necessary for most spacecrafts during initial acquisition. The orbiting 
vehicle leaves the last stage of the booster, and is left with some initial rates at an arbitrary 
3 -axis attitude. The satellite is then stabilized to the sun line, or the local vertical connecting 
the spacecraft's center of gravity and the center of the earth. Only after the initial stabilization 
mode is accomplished does a sensor take over and establish the third positional axis of control. 
This dissertation is concerned with this first mode, i.e., initial acquisition, (or large angle 
control, or restabilization mode); and for purposes of the next chapter the local vertical is as- 
sumed as the reference axis. Thus, interest lies only in the three components cf rate and in the 
three direction cosines, a 13 , a 23 and a 33 , that are formed by the z 3 axis and the spacecraft's 
principal axes. Hence, a five-dimensional problem exists (the three rates &i lt and a> 3 , and 
any two of the three a i3 ), with six state variables involved. The governing equations are 




A 3 - 1 A n 2 

\ 1 2 ) to 3"l + I 2 (“l* "2- " 3 - a l3> a 23' a 3 3 ) 


(hzl^ ^ , 

\ I 3 ) “l“2- + I 3 ("l* <“ 2 . a 13- a 23’ a 3 3 ) 


( 7 . 50 a) 
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(7.50b) 


®13 “3 a 23 “ ^2 a 33 

a 23 " "l a 33 ~ U 3 a 13 

a 3 3 “ "2 a 1 3 " “l a 23 ’ 

which come from Equations (7.45) and (7.47) with the above-stated assumptions. The constraint 
equation used is from Equation (7.24a), 


a J% + a 2 2 3 + a 3% * 1 ■ C 7 ' 51 ) 

Equations (7.50) are in the general form considered in Part I of this dissertation, and only the 
control laws 


Ni = Nj.(S, a 3 ) i = 1,2,3 (7.52) 

need be specified to describe the system completely. 

One additional comment is in order here. Henceforth, cyclical indexed notation is used 
extensively for the sake of brevity. That is, three equations 

x i = y 2 + z 3 

x 2 = y 3 + z i 

X 3 = yj + Z 2 


are written as 


x i “ Vj + z k - 

where the i, 3 , k take on the values of 1,2, 



A graph of the sgn function is shown in Fig- 
ure 7.2 
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Figure 7.2— sgn (x) Function. 



8. INVESTIGATION OF EQUILIBRIUM POINTS 


In this chapter, some characteristic properties of attitude control systems are shown and 
the equilibrium point stability question is investigated. The relationship between potential 
troublesome dynamic behavior and equilibrium point stability is demonstrated. 

One method commonly used for general nonlinear systems' stability investigations, is to 
approximate the system near an equilibrium point linearly, and apply well-defined stability 
criteria. However, for attitude control equations, the linear procedure is usually fruitfully 
pursued for only small regions about the one desired equilibrium state (which of course is 
forced by design to be stable), and at any other points one usually is forced to make invalid as- 
sumptions (with respect to sensors and cross-coupling) in order to obtain practical answers. 
That is, for anything but raw first cuts these methods are unyielding. Perturbation techniques 
cf the type 


x = A x t \f(x) , 


where one tries to extend A = 0 results to the nonlinear case by stretching the parameter A, 
have also not been successful for these highly nonlinear systems that the attitude control en- 
gineer encounters in practice. New techniques are being developed, but at present none can be 
practically utilized. To date, to meet with success one utilizes a computer simulation, and is 
usually forced to make as many Monte Carlo runs as economically feasible. As already pointed 
out, this process is more efficient when the system equilibrium points are all known. 

The equilibrium points for satellite attitude control systems of interest are now determined. 
First the Euler Equations, (7.50a), 

N, = I, OJj (l 2 - I 3 ) eo 2 

N 2 = I 2 «2 - 

# " * « 3 " { 1 l~ I 2) a> l OJ 2 


as developed in the previous chapter, are examined. With R jk = I. - 


Equation (8.1) yields 


N. = I . a>. - Rik wi co, 

ill k 


(8.2) 


where i, j , k = 1, 2, 3 in cyclic order. For the system to be in equilibrium, the rates co i and 
torques N. are time-invariant, thus the rate state variables at equilibrium points are the 
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algebraic solutions to 


N. + 


R jk coj 


( 8 . 3 ) 


or 


That is. 



_ N _ii 




N. R.R.. 

i i ) k i 


( 8 . 4 ) 


This relationship furnishes much information. For any given rigid body, the polarities cf R J2 , 
R 23 , and R 31 are fixed. Hence, the sign combination of N p n„ and N 3 must yield a positive 
quantity for co* . This polarity constraint given by Equation ( 8 . 4 ) immediately reduces the 
number of possible combinations that yield real equilibrium points. (If the control torques de- 
pended upon only the rates, e.g., all sensors are rate gyros, then the three Equations (8.3) or 
( 8 . 4 ), along with the three equations 


N. - IT (<u v cj 2 , a 3 ) , 

completely define the locations of equilibrium states and simplify the problem, since the posi- 
tions do not enter the picture.) 

Further analysis of the a> i components of the equilibrium points is now made. The direction 
cosine Equations (7.50b) which are defined with respect to the local vertical, that is, the a J3 's, 
are 


^13 “ "3 a 23 ~ W 2 a 33 

k 23 = ®33 " "3 a l3 ( 8 - 5 ) 

a33 - ^2 ®13 ~ “l S 23 

There are only two-independent a i3 in (8.5) via the 



i=l 
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constraint. That is, by multiplying a i3 by o> i and adding the three equations, one obtains 


a 13 W 1 + W 2 a 33 ^3 

In the cyclic i , j , k indexed notation, (8.5) is 


a k3 


( 8 . 6 ) 


or 


CO 


k 


a . - + co. a, ,, 
i3 j k3 

a . . 

J 3 


(8.7) 


Substituting (8.7) in the Euler Equation (8.2), 

N i = I i ™i ~ R jk w ■ "k 


(8.2) 


yields 


N. 


I. co. - R., 

ii jk 


a. , + co. a. a 
k3 i j 3 


a i3 +W j a k3 
a j 3 



The above expression, after some algebraic manipulation and use of 



= 1 


reduces to 


N. 

1 

T7 


R jk a f3 a k3 

I . a \ 

i i3 


d . 

R ik dt 


k3 


•ja) 


R., a, a a . a 

jk k3 j3 


21. a. 2 , 

l i3 


( 8 . 8 ) 


Equation (8.8) represents three general equations of motion where the position and velocity 
variables, purposely intertwined, are changing with time. It is shown later that close to an 
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equilibrium point a i3 remains arbitrarily small. (Obviously a i3 = 0 at the equilibrium point.) 
Thus, if the a i3 's are approximately constant, the bracketed terms in (8.8) are small; hence 


N. 

1 

IT 


R jk 3 j 3 a k3 

I . a.\ 
i i3 


(8.9) 


During the time interval for which (8.9) holds, the phase-space trajectories' projection onto the 
three now uncoupled phase planes T. vs. <a. are vertical parabolas. The parabolas are station- 
ary with time if the IT ’s are constant. Which way the parabolas open depends on the sign of 
R jk a j3 a k3 (since E a i3 is positive). The polarity of where the parabolas cross the ordinate 
is the polarity of IT . 

Now several observations are made. The above three uncoupled phase-plane plots repre- 
sent the system's motion if the N. and a i3 are sufficiently constant for a long enough interval 
of time. When these conditions are exactly met, i.e., a i3 - N. =0, then the ok - axis crossings 
of the parabolas yield the values cf that are the equilibrium points cf the entire system (i.e., 
compatible solutions of Equations (7.50) with all derivative terms zero). Thus, for these condi- 
tions it is evident that for positive IT, positive u>. equilibrium points are stable, and negative 
equilibrium points are unstable. These are shown in the first two diagrams of Figure 8.1. For 
negative IT the opposite is true, as shown in the last diagram of the same figure. (Arrows 
indicate direction with time, so that if o>. > 0, then w. increases with time.) That is, if sgn(ir ^ iequib 
= - sgn (R jk a. 3 a k3 ^ iequib ) = +1 (- 1), then the equilibrium point is stable (unstable). (Eater on, 
a torque bang-bang system is discussed in which two unstable equilibrium points trap a phase 
trajectory between them, and around a third point, such that a limit cycle is generated.) Even 
though each parabola may contain two equilibrium points, one cf these points is usually elimin- 
ated by the control laws. F a parabola opens up (down) and N. is positive (negative), then no 
equilibrium points exist, and the algebraic system will yield only non-real solutions. This 
occurs when sgn(R„ a j3 a k3 N.)= +1. Thus, sgn (R jk a j3 a k3 it) = -1 is a necessary condition 
f° r "iequib to exist. If N. is zero, then there exists only the one quasi- stable equilibrium point. 
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(Note: Many more polarity constraints can be derived for a i3 , and N. to help locate equi- 
librium points. These are treated later, and tabulated in Table 8.2.) 

The assumption that a i3 remains sufficiently small around an equilibrium point is now in- 
vestigated. The three equations 


a 


13 


a j3 


co. a 


k3 


| i, j, k, - 1,2, 3 cyclic indices 


( 8 . 10 ) 


lead to a single uncoupled third-order equation for each a. . As previously noted, the above 
three equations are not independent; the algebraic relationship 



( 8 . 11 ) 


must be satisfied. After multiplying each a i3 in Equations (8.10 ) by co i and adding, one obtains 


®13 "l + *23 U 2 + a 33 "3 = 0 - (8.12) 

Taking the second derivative of Equations (8.10) with respect to time yields, after some algebraic 
manipulation, 


+ 


(co. 2 + co 2 + co 2 ) k i3 - - I a i3 ^ (co 2 + co 2 ) 



r d ,. 

. a 


r d /• \ • i 

+ a i3 

Ldt K + "i “v 

+ "x "j. 

a k3 

Ldt (“j ~ oj i "k) _ "i "kJ 


(8.13) 


Like Equation (8.8), the above equation represents purposely intertwined general equations of 
motion in which the state variables are changing with time. However, the equations are un- 
coupled when the cods are assumed constant (which they are at an equilibrium point); in this 
case, the entire right-hand side of Equations (8.13) is zero. That is, 

a i3 1 " 2 a i3 ~ 0 » i = 1, 2, 3 (8.14) 

where 


co - y co 2 + co 2 + co 2 - constant (8.15) 

is the magnitude of the system's angular velocity vector. The three simple third-order linear 
Equations (8.14) are solved to yield the behavior cf a i3 when w. = co. = co k =0. That is, 

a J3 (t) = A 3 + B. sinwt + Cf cos wt . (8.16) 
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Here, the A i5 B i 5 and C i are constants, but not arbitrary constants. They are determined by 
the initial conditions and the constraint (8.11). The relationship cf the nine constants is in- 
vestigated in the Appendix. The important point here is that by letting all acceleration terms 
vanish, the three a i3 's are decoupled and can now each be examined independently, using, as 
an aid, phase plane plots. In each phase plane, 

a. 2 

( a i3 _A i ) 2 + ^F = constant = K 2 , (817) 


and the phase portrait consists of a family of periodic trajectories that are concentric ellipses 
centered at the equilibrium point, a 13 = Aj. Hence, it is clear from (8.17) and Figure 8.2, that 

starting arbitrarily close to the centroid at 
A., the trajectory hovers around it, with the 
maximum a i3 = Kw occuring at a. 3 = A i . 

Thus, for sufficiently small K, a i3 stays ar- 
bitrarily small. 

The above discussion is now summarized. 
At an equilibrium point, all dynamic quanti- 
3 ties remain constant with respect to time. 
When it is assumed that near an equilibrium 
point the a i3 = 0 but not the oi i , then, de- 
pending on certain polarity considerations, 

Figure 8.2- Direction Cosine Phase Plane one concludes that the state cf the system 

Portrait for A = 0. moves either toward or away from the equi- 

librium state. When it is assumed that the 

A = 0 but not the a i3 , and when the system is started near an a i3 = 0 system equilibrium 
state, then the system state remains near the equilibrium state. Hence, if the system state is 
very close to an equilibrium point at t =0, then it is the A-vs.-gj. parabola polarity that de- 
termines stability in some neighborhood of the equilibrium point. This conclusion is fortified 
by the fact that usually most main torque generators used in attitude control systems are used 
in an on-off fashion. Therefore, unless a torquer is duty-cycling (as discussed in the next 
chapter), its bang-bang characteristic gives it a digitized constant level in a neighborhood of 
cf sensor signals such that N. is constant, thereby keeping the A-vs.-tu. parabolas almost 

stationary with time. For example, an ideal 




pneumatic loop, consisting of voltage from 
a rate gyroscope feeding into a double-ended 
threshold detector that drives pneumatic 
thruster solenoids, is shown in block diagram 
form in Figure 8.3. 
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Figure 8.3— Ideal Pneumatic Loop. 


From Equation (8.4), it is seen that if 
the N. 's take on only discrete levels, then 




out”the system so that it correlates with the 
mathematical system described in Part I of 
this dissertation, the bang-off-bang pneu- 
matic controller is roundedout to that shown 
in Figure 8.4. This continuous characteristic 
is: 



Further strengthening cf the conclusion that the values of oj. (and the control laws which 
yield Nj) at an equilibrium point completely determine the stability of the equilibrium point is 
obtained by an approximate linear expansion of the equations of motion about that equilibrium 
point. Defining the values of the state variables at the equilibrium point as <u Ie , a 2et co 3e , a 13e , 
a 23e , a 33e , the transformation, 


"i = u i + "le 

"2 = U 2 + W 2e 

W 3 = U 3 + U 3e 

(8.18) 

a !3 _ U 4 + a i3e 

a 23 = U 5 + 3 23e 

a 33 = U 6 + a 33e 

converts Equations (8.1) and (8,5), for w near equilibrium, into 

N i % I iV( I r I k)he u » t \V B iA) ( 8 - 19 > 

for cyclic i, j, k = 1, 2, 3, and 

= ( U k + "ke) ( U m + a j3e) " ( U j +CJ je) ( U n + a k3e) 

~ "ke U m + a j3e U k + "ke al3e " &k3e ^ ^ ^ ^ 

for cyclic l, m, n = 4, 5, 6, where all higher order nonlinear terms in u are neglected. 
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Equations (8.19) and (8.20) are equivalent to the linear matrix equation 


r , ~\ 

U, 




33e 

_a 23e 


“33 e 

0 

S 13e 


“23e 


13e 


0 

0 

0 

0 

-<y, 

3e 


0 

0 

0 

“ 3 e 

0 


0 

0 

0 

™2e 

“le 

0 


r ' 
u, 


r. 


N. 


1 + 


N„ 


+ <T 3l“3e“le 


+ Cr 12"le W 2e 


“le a 33e “3e 3 13e 


rs.2i; 


-V 


where a. ti = (i. -I.yi k . When the matrices are partitioned as indicated in (8.21), the roles cf 
Equations (7. 37) through (7.39), (7. 44) and (7.47) are clear and (8.21) is derivable from these 
equations with the proper assumptions. The linear Equation (8.21) has the form 


(u) = [D] (u) + (E) , (8.22) 

where (u) is the variable-column matrix, and [D] and (E) are constants. It is important to 
note that 


det [Dl = |d| = 0 . ( 8 . 23 ) 

From a stability viewpoint, if the eigenvalue problem is attacked by solving 

| [D] -A. [U] | =0, (8.24) 

where [U] is the six-by-six identity matrix, the following sixth-order polynomial in A is 
obtained: 



~k 

^23 “3e 

^23 “2e 

M“le + “ 2 e +C0 3l +A - 2 ) 

a 31 “ 3 e 

-A 

^31 “le 


ff 12 “ 2 e 

a X2 “le 

-k 


0 , 


(8.25) 


or 


^ + He +0) 2e +«£)] 


12 


,2 

'3l“le 


J 12 (J 23° J 2e 




, CO? CO? CO ? 
le 2e 3e 


(8.26) 
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Several observations are now made concerning the characteristic Equation (8.26). Firstly, 
nowhere do the direction cosines a i3e appear. That is, all six eigenvalues (X-roots) are inde- 
pendent of the spacecraft attitude near the equilibrium point. (A constant value cf N ; is of 
course assumed. That is, N. is not a function cf u (or time) and hence does not enter into [D] 
but only into (E) in Equation (8.22). Thus, by letting (v) = (u) + [D] -1 (E), Equation (8.22) be- 
comes (+) = [D] ( v) and the eigenvalues are the same.) Another interesting observation is that 
three cf the characteristic roots of Equation (8.26) are 


k l ~ 0 ' k 2 “ ^ > 




where 


CO 


e 



We 


+ CO 


2 

3e 


is a real number (and the same frequency as in Equation (8.16) where 


CO. - CO. 


was held con- 


-\[U] 


= 0, where J is as in 


stant). These are the three eigenvalues obtained from 
Equation (7.39). Thus, the linear analysis does not yield the complete stability picture and the 
nonlinear terms must be accounted for. It is also interesting to note that the remaining three 
roots, \ 4 , k sf and \ 6 , depend on the spacecraft moments cf inertia as well as the co ie ’s for their 
realness as well as for their polarity. 


Another aspect of the attitude control system's equilibrium points is that cf polarity com- 
binations mentioned earlier in this chapter. For the sake cf definitiveness, it is assumed that 
the spacecraft has the following inertia relationships: 


I, > I 2 > I, > 0 . ( 8 .27) 

(Noloss in generality results from this assumption since similar results are obtained for any 
inertia magnitude order.) Hence, 


R 23 > 0 (8.28) 

*31 < 0 


and Equations (8.3) become 

Nj = - Kl 012 0)3 

N 2 = K 2 Wj o> 3 (8.29) 

N 3 ~ _ K 3 "l "2 ' 
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where K 3 (i = 1 , 2 , 3 ) is positive. The above three equations, plus Equations ( 8 . 5 ) with no 
derivative terms, i.e., 


23 

“2 a33 


33 

= "3®13 

( 8 . 30 ) 

13 

= u t a 23 



are satisfied at any equilibrium point. In addition, 


* 1 2 3 + + * 3 \ = 1 < 8 - 31 ) 

must be true. There are eight possible polarity combinations for the three &>., and from Equa- 
tions ( 8 . 29 ) each of these combinations yields a specific polarity for N. . They are shown in 
Table 8 . 1 . 


Table 8.1 


Combination 

"1 

"2 

— 

“3 

N i 

N 2 

N 3 

1 

+ 

+ 

+ 

- 

+ 

— 

2 

+ 


- 

+ 

- 

- 

3 

+ 

- 

+ 

+ 

+ 

+ 

4 

+ 

- 

- 

- 

- 

+ 

5 

- 

+ 


- 

- 

+ 

6 

- 

+ 

- 

+ 

+ 

+ 

7 

- 

- 

+ 


- 

- 

8 

- 

- 

- 

- 

+ 

- 


Table 8.2 



From Equations (8.4), 

S gn(N 1 N 2 N 3 ) - +1 ( 8 . 32 ) 

for ccj. 2 positive (and thus &j. real). (Note 
that 

sgn(N j N k /N i ) - sgn(N lN2 N 3 ) , 

and 

sgn (R-k/Rij R ki ) = -1 

from Equation (8.28).) Equation ( 8 . 32 ) is, of 
course, true for all combinations in Table 8.1 
In total, there are 64 polarity combinations 
cf all the variables, that is eight direction 
cosine (a 13 , a 23 , a 33 ) combinations for every 
combination of rates in Table 8 . 1 . However, 
Equations ( 8 . 30 ) reduce the number of cases 
to 16 as given in Table 8 . 2 . 

Table 8.2 shows all possible combina- 
tions of parameters for arbitrary control 
laws (assu’ming a spacecraft with inertias as 
in Equation ( 8 . 27 )) for any real equilibrium 
points. It is evident from the table that there 
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are only four permissible control torque po- 
larity combinations. Each one exists in only 
two octants that are mirror images* of each 
other on the direction cosine unity sphere. 
The unity sphere is a useful geometric de- 
vice for displaying attitude, and is defined as 
a sphere in a 13 , a 23 , a 33 space cf unit radius 
and center at the origin. (See Figure 8.5). 
Allattitudes of the spacecraft lie on the sur- 
face cf this sphere since 


D 

i= 1 


:. 2 , = 
i3 



Each surface quadrant cf the direction cosine 
unity sphere allows equilibrium points to 
exist only in two (mirror image) octants of 
three-dimensional co 1} co 2 , a> 3 rate space. The 

two mirror image points fall out from the two roots of Equation ( 8 . 4 ). Table 8.2 and its inter- 
pretation (which comes with practice) is an extremely useful analytical tool in locating equi- 
librium points in a system. 


The above results are independent cf any specific attitude control system with explicit con- 
trol laws. This chapter has been an attempt at understanding equilibrium points that exist in 
an attitude control system. This analysis, coupled with the mathematical development cf Part I 
that puts equilibrium points inside limit cycles, establishes a basis for finding potentially 
troublesome regions of the system. In the following chapter, examples cf practical application 
are given in the hope cf demonstrating the power behind this type cf an approach when utilized 
with a computer simulation. 


*Miiror image points in 3-dimensional space is where the polarity of each component of one point is the negative of its respective 
component of the other point, i.e., the position vector to both points are equal and opposite. 
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9. PRACTICAL APPLICATIONS 


The equilibrium analysis technique de- 
veloped in the previous chapters has been 
successfully applied to actual complex atti- 
tude control systems. The power of this ap- 
proach is demonstrated by the resulting 
usefulness from these practical applications. 

Thus, in this chapter, the method is illus- 
trated by considering the following systems: 

I) The NASA* advanced meteorological 
satellite. Nimbus I, (Figure 9.1), an earth- 
oriented spacecraft. The analysis shows the 
shortcomings of the attitude control system 
for large perturbations (such as those which 
wouldnecessitate reacquisition of the earth's 
local vertical) or large initial conditions at 
the time of separation from the launch ve- 
hicle. Nimbus I was launched August 28 , 

1964 , and had a fairly successful active life 
until September 23 , 1964 . Its solar array drive mechanism then locked up, preventing efficient 
energy conversion from the sun. As a result, the spacecraft experienced intermittent power 
outages. This, in turn, led to erroneous signals in the pneumatic torquing system. The er- 
roneous error signals caused the vehicle to spin at such a high angular velocity that even if the 
power situation had corrected itself, the spacecraft could not be restabilized (Reference 26 ) . 

ID The follow-on Nimbus weather satellites II and 111 which have progressive modifica- 
tions in the attitude control systems. Although improved, due to cost limitations involved in 
change, they still have similar stability-in-the-large potential dangers. Nimbus llwas launched 
May 15, 1966 . The initial, conditions for earth acquisition were small and the spacecraft stabil- 
ized successfully. The control system is still performing satisfactorily. Nimbus lllwill be 
launched in the late part of 1967 . It is interesting to note that several easily implemented 
schemes were proposed for remedying certain non-stable characteristics of the system and at 
first, after several hundred thousand computer runs, these schemes appeared attractive. How- 
ever, after a thorough equilibrium analysis (Reference 27), it was shown that in several cases 
the original equilibrium points are moved to different locations in state space, bringing their 
unstable behavior with them. Also, new equilibrium points are introduced. The large number 



^National Aeronautics & Space Administration. 
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of computer runs never located these isolated trouble areas. As an example, runs are made 
where the initial attitude and the roll and pitch rates are held constant, while the yaw initial 
rate is varied as a parameter for each simulated run. Later, when equilibrium points are 
calculated, one is found to lie in between two of these yaw-rate runs. When the computer is 
initialized at the equilibrium point, the system is unstable; it "meanders" so to speak (in al- 
most periodic fashion), ejecting gas, and then finally takes off. Here is a "hole" in a region of 
state space declared a stable region because cf the results from a grid of computer runs in the 
region. This example illustrates that unless the initial condition grid is infinitesimally small 
(whichtakes an infinitely large number of runs, especially for high-dimensional systems), a 
brute-force coverage pattern or a Monte-Carlo random approach does not yield a high degree 
of confidence that the computer will find all the "holes" of a complex nonlinear system. The 
equilibrium analysis is designed to substantially boost this confidence factor. 

Ill) The advanced Nimbus satellite (Nimbus D) which has a completely new set of attitude 
control laws and sensors. This vehicle is scheduled for launch in 1969. The equilibrium point 
approach was used in the original design (which is now being completed) to eliminate all equi- 
librium points in a reasonably bounded region of interest of state space. This application is 
described in more detail later in this chapter. 

The equilibrium approach is to be attempted for other major scientific NASA spacecraft, 
such as the Orbiting Astronomical Observatory (OAO), which is an inertially-oriented satellite 
to study the stars; the Orbiting Geophysical Observatory (OGO), which is earth-oriented and 

studies many physical phenomena in a region 
of space around the earth; the Aerobee Sound- 
ing Rocket 350, which is a short-lived eco- 
nomical inertially-oriented spacecraft that is 
flexible so as to study many items of inter- 
est;the Applied Technology Satellite (ATS-F) 
whose design has just commenced and will be 
used to demonstrate space technologyin many 
areas; and possibly many others (Orbiting 
Solar Observatory, Deep Space and Galactic 
Probes, Apollo Application Missions, etc.). 

A brief background description of the 
Nimbus weather satellite is now presented. 
How the equilibrium approach is utilized, 
and the results, are given. (For more de- 
tailed description of Nimbus, see Reference 
28, 29.) Nimbus I is shown in Figure 9.2. 
The spacecraft stands 9-1/2 feet tall, weighs 
850 pounds, and has moments of inertia 
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Figure 9.2 — Nimbus I. 
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approximately equal to: 


I, = 222 slug-ft 2 * , I 2 - 19£>slug-ft 2 , I, = 90 slug-ft 2 , (9-1) 

and peak-power capability of approximately 470 watts generated from a solar-cell array. (It is 
noted that many simplifying assumptions are made for the sake of brevity. However, the actual 
detailed analysis includes the entire physical system as well as it can be determined; all as- 
sumptions are looked at qualitatively if not quantitatively to see their effects on the final re- 
sults. One such assumption, aside from taking the structure as a rigid body, is to ignore the 
smalltime-varying property of the moments of inertia due to the large solar- array paddles ro- 
tating to follow the sun.) The spacecraft is broken down into three major physical sections: 

1) The sensory ring on the bottom which houses all the meteorological equipment that looks 
down at the earth, (such as T.V. cameras, infrared energy sensors for obtaining three- 
dimensional pictures of clouds, etc.), and the associated subsystems such as batteries and power 
supply, communications and antennas, spacecraft clock, etc. 

2) Two 3-1/4 x 8 feet solar array pad- 
'dles that are driven via a closed control loop 
to follow the sun for maximum energy input 
(Figure 9.3). The orbit is a nominal 500- 
nautical-mile-altitude, circular, "high-noon," 
nearly polar orbit. That is, launch time is 
local midnight, and the inclination angle of 
the orbit (^98.7 degrees to the equator) is 
such that the orbital plane precesses one 
revolution per year, at the same rate as the 
earth moving around the sun f (Figure 9.5). 

Thus, the sun remains nominally in the or- 
bital plane (i.e., the orbital plane is sun- 
synchronous), and, if the attitude control sys- 
tem aligns the vehicle in the orbital plane, 
then the solar arraydrive needs only one de- 
gree of freedom (about the orbital pitch axis). 

The solar array drive control system is a 
straightforward feedback loop and is func- 
tionally designed and analyzed in a simple 
classical manner, but the mechanical hard- 
ware problems, as demonstrated by the first 
Nimbus flight failure, are not so simple. 



*3 



Figure 9.3 — Nimbus Spacecraft Body Axes. 


*1 slug-ft 2 = 1 ft-lb-sec 2 . 

tThis orbit also affords the spacecraft cameras complete earth coverage twice a day, during local noon and midnight, as shown in 
Figure 9.4. 




Figure 9.4 — Nimbus Picture Coverage. 
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Figure 9.5 — Orbit-Plane Regression (Autumnal Equinox). 


3) The attitude control system and all 
its associated subsystems which, of course, 
are the parts of the satellite of most con- 
cern for this dissertation. After a Thor- 
Agena combination boosts and injects Nim- 
bus into orbit, the satellite separates from 
the Agena at approximately a vertical po- 
sition with angular rates (hopefully) less 
than ldegree/sec in all three axes. The 
solar paddles are unfolded, and the control 
system is initiated. It brings the spacecraft 
from its initial conditions to within a small 
neighborhood cf the desired equilibrium 
point in a short time, and with a small per- 
centage of available control gas. The de- 
sired equilibrium point expressed in the 
notation cf Chapter 7 is: 
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(9.2) 


where (X) = [A] (y) and (y) = [a] (z). ff the acquisition and control response times are small 
(on the order of seconds) compared to the orbital period (r = 2 rr/Q 0 for a 500 nautical mile al- 
titude -103.5 minutes) as they are in this case, orbital rate may be neglected compared to at- 
titude rates, so that 

M = to] , (9.3) 

and the system is as described in Chapter 8. Equation (9.2) then describes a satellite whose 
vertical (yaw) axis is aligned to the earth's local vertical as the satellite orbits round the 
earth, and whose body-roll and yaw axes lie in the orbital plane. That is, the satellite is in- 
ertially rotating at one revolution per orbit about an axis normal to the orbital plane. The ac- 
curacy and rate specifications for the control system are that the angular errors (in an Euler 
small-angle sense) be less than one degree, and the angular error rates be less than 0.05 de- 
grees per second with respect to the orbital reference axes. 

To meet these requirements, Nimbus I has as control system sensors: 

A) Two infrared horizon scanners which ideally provide sufficient information to determine 
the location of the earth's local vertical. 

B) A rate gyroscope used as a gyro-compass to furnish yaw control. 

C) Sun sensors for backup yaw control to the sun line. 

For control actuators (muscles), there are: 

D) A pneumatic cold gas mass expulsion system capable of changing the system momentum. 

E) Reaction flywheels, which are two-phase linearly wound inside-out motors (heavy pan- 
cake rotor revolves around inside stator), capable cf storing momentum. 

Before the control laws that connect the sensors and muscles are defined, a brief descrip- 
tion is stated for each. One must realize that volumes could be (and have been) written (e.g. 
References 30, 31) on any one particular area of an attitude control system, such as a specific 
scanner, but here the intent is to oversimplify for background familiarity. First the sensors 
are discussed. 


Sensors: 

A) Each horizon scanner has an optical prism that is motor-driven to sweep out a conical 
field of view, with the cone axis along the body-roll axis (x 1 ). There is one scanner fore and 
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FIELD OF VIEW SCANNER ROTATION 



Figure 9.6a— Earth as Seen by Scanners in 
Stabilized Nimbus. 
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Figure 9.6b— Roll and Pitch Error Computation. 


another aft of the spacecraft (Figure 9.6a). 
Incident infrared (IR) energy is detected by 
an immersed thermistor bolometer, and a 
voltage pulse is obtained as the scan path 
moves through cold space to hot earth and 
back to cold space again. If the spacecraft 
has a roll error, then a vertical reference 
pulse generated at the bottom of the scan 
cone (fixed with respect to the spacecraft), 
is displaced from the center of the earth IR 
pulse, and this displacement is a measure of 
the roll error. To obtain pitch information, 
a scanner cannot be used along the pitch 
axis, since the solar paddles block the field 
of view; therefore, two roll scanners are 
used back-to-back. A pitch error angle 
means that one scanner apparently sees a 
larger earth than the other. Hence, the dif- 
ference in the two IR earth pulse widths is a 
measure of the pitch error (Figure 9.6b). 

The scanners are extremely nonlinear roll- 
pitch cross-coupled error detectors for all 
but very small angles. This is due to (l)the 
geometry of the situation— for instance, if 
the spacecraft pitches far enough one scan- 
ner leaves the earth. Pitching still further 
until the spacecraft is almost horizontal, the 
other scanner cone sees only earth; (2) the 
digital error processing must be severely 
constrained— for example, not only are sun- 
shades used, but the bolometer output signal 
is electrically clamped to zero at the top of 
the scan cone; this ensures that if the sun 
comes into the field of view (the sun is a 
very hot IR source) it is not processed as is 
the earth (which is actually down below); (3) 
the earth is not an ideal homogeneous spher- 
ical IR source — for example, high cold clouds 
next to a warm patch on earth may lead the 
sensor to suppose that a cold-space/hot-earth 
horizon is falsely within the earth disc, thus 
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producing significant sensing errors; (4) the scanners and processors are non-ideal pieces of 
hardware in many physical respects. Thus, it is a formidable task to mathematically repre- 
sent realistic roll and pitch control errors as 


^1 ^1 ( a i3’ a 23 ’ 3 3 3 ) 

^2 ~ V 2 ( a !3 ’ a 23’ a 33) 


(9.4) 


and one is compelled to use a graphical representation (Reference 32 , also page 9.11 and Fig- 
ures 9. 10through9. 13 later in this chapter). The graphs consist of a family cf parametric curves 
on the a i3 (i = 1 , 2 , 3) direction-cosine unity sphere that are projected on the a 13 , a 23 plane 
for positive and negative a 33 . A highly complex hardware simulation that includes the actual 
geometry can also be used statically to find the errors needed in the equilibrium point search. 
Note that for an assumed spherical earth, V, and V 2 do not depend on yaw, and hence are only 
functions of a i3 (i = 1 , 2 , 3 ). 

B) To understand the gyro-compass error detector, assume that the spacecraft is per- 
fectly aligned with the orbital reference axes with zero relative rate. Then, the only rate a 
gyro measures is the inertial rate, perpendicular to the orbital plane, of the reference axes 
going around the earth, i.e., -fi 0 y 2 . Hence, if the input axis of the gyro is along the body-roll 
axis x 1? any misalignment of x ( out of the orbital plane causes the gyro to sense a component 
of orbital rate proportional to the sine cf the misalignment yaw angle. Thus, the gyro rate out- 
put is a measure of the yaw error and is used as such. The input axis of the gyro is tilted up 
toward the -x 3 axis by a small angle y so that the gyro output voltage also contains a small 
component cf yaw error rate, « 3 , for damping cf the yaw control loop. The yaw error voltage 
is 


V, - Wj cosy - W 3 siny (9.5) 

in the notation of Chapter 7. If = 0 is assumed, then V, no longer contains yaw angle in- 
formation, but only undesired roll rate and desired yaw rate. That is. 


V, - a> l cosy - co 3 siny , 


(9.6) 


and the system exhibits no yaw positional control. The gyro-compass operation is illustrated 
in Figure 9.7 for zero roll and pitch errors. 

C) The sun sensor is used in a backup failure mode and is not considered in this chapter. 


88 




g; 

w = - n 0 ? 2 fc 'Ay 3 = ■ ( Sin '/’ J 1 + co s 'A’ ? 2 ) + ^*3 


Muscles: 

D) In Nimbus, bang-off-bang controllers 
for the pneumatics are used. Current from 
a power stage of a double-ended threshold de- 
tector with hysteresis is fed to the proper 
gas solenoid valve which, when open, pro- 
duces a thrust. This thrust is collinear with 
the spacecraft center of mass; hence, an ex- 
ternal torque is produced. A simplified sche- 
matic representation is shown in Figure 
9.8(a). In Figure 9.8(b), the system is further 
simplified by ignoring hysteresis, valve trans- 
port delays, pressure buildup lags, etc., and 
combining everything in a smooth pseudo- 
bang-off-bang curve with sufficiently sharp 
slopes to be realistic enough for present pur- 
poses. One of these torque generators exists 
for each of the three axes: roll, pitch, and 
yaw. 




- 0 . sin \p , W, fi cos \p , W 3 - 


E) Flywheels are used primarily for 
fine pointing control for small and cyclic- 
type disturbances. (For cumulative disturb- 
ances, the flywheels saturate and become 
ineffectual; thus, gas is fired to unload mo- 
mentum when the wheels reach a certain speed). The wheel torque levels are much smaller than 
the gas system; therefore the wheels are neglected for present purposes. 


thus, V 3 — — (fl 0 cos y ) sin </r - (sin , y)i/f and, if 

fl Q = 0 (i.e. % = 5 j), then V, contains no yaw position information 

Figure 9.7— Gyro-Compass Operation. 


The system is connected together by electronic implementation of control laws, involving 
loop compensation networks designed to give the loops proper response. When searching for 
equilibrium points, one looks for time-invariant behavior; therefore noise filters, lead-lag 
networks, etc., need not be considered dynamically. Only the DC steady-state gains are at 
issue. Limiting, and physical saturation are also considered. The control laws of Nimbus I re- 
duce to those statically shown in Figure 9.9. However, early in the design it was realized, 
and vividly dramatized by computer runs, that the system does not work well for all reasonable 
rates. The reason is the yaw loop responds to five times more roll rate than yaw rate, so that 
the yaw gas produces still higher yaw rates. These rates couple back through the dynamics as 
roll and pitch torques, as evidenced by Equation 8.1. The loop is designed for fine performance 
as well as coarse; so the V 3 control law is needed to obtain yaw-angle information which is 
only contained in the oq term.) The yaw gas is throttled down by allowing it to pulse for 
at most one half second out cf every thirty seconds. Thus, N 3 (compared to N, andN 2 ) is 
small for the first few important seconds of action. Hence, for present purposes, N 3 = 0 
is assumed. 
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+ ROLL THRUST 


-ROLL THRUST 


TORQUE Ni 



|Ni|« Tj sgn (| vj | -M-,) 


Figure 9.8— Simpl ified Pneumatic System. 

The work is now made simple by the groundwork already laid by this and previous chapters. 
All equilibrium points are solutions to the following set of relationships: 


^ K 1 ( V 1 ( a 13> a 23> a 33))) ( I 2 I3) a) 2 C °3 

(9.7a-l) 

^K 2 (V 2 ( a 13- a 23> a 33))) = "(is -I l) " 3 “l 

(9. 7a- 2) 

= 0 = '( I 1 _I 2 ) “i“i 

(9.„7a-3) 

W 3 a 23 = "2 a 33 

(9.7b-l) 

“l a 3 3 ~ W 3 a 13 

(9.7b-2) 

W 2 a 1 3 ~ "l 3 2 3 

(9.7b-3) 

a 13 + a 23 + a 3 2 3 ~ 1 ' 

(9.7c) 
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TORQUE ERROR SIGNALS, v, = K, (Vi), i= 1, 2, 3, 
v 1M = v 2M = 12“ SATURATION FOR NIMBUS I 


Figure 9.9-Simpl ified Static Control Laws 


V, (a,,, a 23 > a 33) an d v - < a ”- a 23 > a 33) are gi yen graphically in Figure 9.10 through 9 . 13 . 

(Vj (0, 0, 1) = V 2 (0, 0, 1) = Oby design, to make the origin, - u 2 = - a 13 - a 23 - 0 , 

«=a 33 : 1 an equilibrium point.) On these curves, v 2 = K, (Vj)’ and v 2 = K, (V 2 ) (which include 
all saturations) are indicated by crossed-hatched regions. Thus, N, and N, are defined. Since 
a 33 is redundant via Equation ( 9 . 7 c), the curves are presented as two projections within the 
a j 3 -versus- a„-plane unity circles, one for positive a 33 and one for negative a 33 . Therefore, 
only four sets of curves in total are necessary. How they are generated (Reference 32) is in- 
cidental to the discussion. The curves are used as follows: for any (a,,, a 23 ), the family of 
curves (interpolation is needed between curves) have numbers which are in units of degrees, 
representing the signals fed to the gas bang-off-bang controllers, i.e., v. = Kj (a,,, a 23 , 

a 3 5)) with’the crossed-hatched regions representing saturated signals. There are shaded 
clamped and indeterminant regions that, for example, represent a scan cone completely on the 
earth so that the error signal depends on a random cold-cloud distribution over the earth at 
that particular time. Hence, in these regions, the signals are anywhere from minus saturation 
to plus saturation. 
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/ ////\ 12° SATURATION REGION 
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NUMBERS WITHOUT DECIMALS 

ARE Vi (a 13 , ° 23 / ^33) 
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Figure 9.10— Roll Error Versus Attitude, a 33 > 0. 
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12° SATURATION REGIONS ARE NOT SHOWN 
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12' SATURATION REGIONS ARE NOT SHOWN 
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IN DEGREES 
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Equations (9.7) are now solved. From (9.7a-3), either or o> 2 (or both) are zero. If cu 1 = o, 
then N 2 = o (therefore |k 2 (V 2 )j <M 2 ) from (9.7a-2). From Equations (9.7b) either cq, = aq - 0 
(therefore N, = 0), or ce 2 = a 13 = 0 (therefore N, - 0), or co 3 - a 13 = 0 (therefore N, = 0), or only 
a !3 is zero. The various cases are traced out in Figure 9.14. There are only six distinct 
categories, as tabulated in Figure 9.15. 

Case I, which includes the origin, indicates that if all the rates are zero, any attitude for 
which the scanners have near-zero output are equilibrium points. From the K. (V;) curves, 
there are many regions where this is true. The only. way to eliminate these cases, is to change 
the sensors and make the deadbands zero. Cases IF 111 and IV allow co L to have any value when 


n 3 =o 


Oi = N 2 = 0 


ij 2 u 3 0 u 2 — Q13— 0 u 3 — a^ 3 — 0 a 13 -0 

Nj,=0 Nj = 0 1^=0 
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Figure 9.14— Equilibrium Cases. 
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Figure 9.15— Table of Equilibrium Cases. 
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a i3 = +i (i = 1, 2, 3), as long as the sensors have near- zero outputs. Therefore, Case 11 ad- 
mits a right-side-up or upside-down spacecraft which is spinning about the yaw (x 3 ) body axis. 
Case 111 cannot exist at a 23 = ti, since K, (v, (0, + 1 , 0))<~Mj. However, owing to nonsymmet- 
rical scanners, the a 23 = -l point can exist. In caseTV, the vehicle is pitched 90 degrees and 
rolling at an arbitrary rate about its plus or minus roll axis which is vertical. 

Cases I through IV arise from false nulls of the sensors, that is, near- zero outputs at at- 
titudes other than the origin. This is a sensor problem that is not easy to eliminate. However, 
in all four cases, no gas is used (i.e., N. = 0, i = 1,2, 3). Hence, there is the comforting 
fact that the spacecraft is not lost (i.e., gas is not expended) while the long-term effects (gyro- 
compass, flywheels and other small effects which have been neglected) pull the system away 
from these points even if they turn out to be stable. These are easy points to put into the full- 
blown computer simulation, and the behavior can be observed. In fact, these are the obvious 
points associated with any physical system with non-ideal sensors, and the points are usually 
easily found. If these points produce instabilities, then there exists the constraint of non-all- 
attitude stabilization, even from zero rates. Such turns out to be the case for Nimbus I. 

It is the more subtle and interesting cases V and VI that are easily overlooked. Here the 
sensors do not have to have zero outputs; yet the vehicle is in equilibrium and using up gas. 
These are the "nasty" points that cause trouble in a system design. For illustrative purposes, 
only Case V is examined, since Case VI is analogous. First, K, (V 2 ) curves are used, and only 
the regions of a 13 , a 23 where the j K 2 ( V 2 ) | curves are <M 2 are considered. These regions are 
the a 23 axes and the three shaded areas shown in Figures 9.12 and 9.13. In these regions, if 
o> 1 = 0, which it does as seen from Equation (8.4), then Equations (9.7a-2 and -3) are satisfied. 
But for Case V, a 13 = 0 is necessary to satisfy Equations (9.7b-2 and -3) as indicated in Fig- 
ure 9.15. Thus in Figures 9.12 and 9.13 only the a 23 axes need be considered, The remaining 
relationships from Equations (9.7) are 


2 - 
a 3 3 

1 - 

a 2 
a 23 

(9.8a) 

* *3 a 23 

— CO 

2 a 33 

(9.8b) 

Ni (k, (v, (o, a 23 , ±yr 

-a, 2 ,' 

23; 

~ ~ (*2 ~ * 3 ) W 2 ^3 ‘ 

(9.8c) 

The control laws approximately yield 




Nj = 0 

for 

|Ki M < Mj 


N, = T, 

for 

K, (V,) > Mj 


N, = -Tj 

for 

K, (V 3 ) < “Mj 

(9.9) 

0 < N, < Tj 

for 

K, (V,) = M, 


-Ti < N, < 0 

for 

K, (V,) = ~Mj , 
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as shown in Figure 9.9. Now an entire set of equilibrium points is obtained. First, plotting 
K, (V,) from the K, (vj curves with a 13 = 0, and then applying Equation (9.9), yields the 
N, = N, (a 23 ) shown in Figure 9.16. Clearly, 

sgn(N 1 a 23 ) = -1 ( 9 . 10 ) 

for all but the indeterminant region. Combining Equations (9.8) yields 

( a 2 3 ) a 2 3 /®23 a 33 ' (9-11) 


where R 23 = I, - I 3 > 0 for Nimbus. Therefore, aside from the indeterminant region, real values 
of a i 2 exist when a 33 > 0, that is, the vehicle is right-side-up. In the indeterminant region, 
sgn(a 23 a 33 ) = + 1; hence, N, < 0 for an equilibrium to exist. Therefore, a large set of equilib- 
rium points are found in the deterministic right-side-up region. They are: 




a 23 - arbitrary 


*33 = + yi ~ a 2 2 3 


oij = 0 

= "2 ( a 2 3 ) ^ 9 ‘ 12) 
W 3 = ( a 3 3 / a 2 3 ) "2 

N, = Ma 23 ) 

N, = 0 

N 3 = 0 , 

'where N x (a 23 )is given in Figure 9.16 and 



( 9 . 13 ) 


which is plotted in Figure 9.17 with Nimbus Inumbers. Physically interpreting this equilibrium 
point, one concludes for a 13 = 0 that the body-roll axis remains perpendicular to the local 
vertical y 3 . Hence, as a 23 varies while a 13 = 0, the spacecraft is rolled about the body-roll 
axis. In fact, for the yaw, roll, pitch Euler angle sequence such as described in Chapter 6, the 
a 23 direction cosine is expressed as 


a 23 = s in 4> , 


(9.14) 
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where <t> is the roll angle. Thus, rolling the 
spacecraft by an angle 4> generates a signal 
to fire roll gas, and the roll torque isjust 
enough to counteract the pitch-yaw coupling 
torque so that co t remains zero and all state 
variables are constant. 

The stability of this equilibrium point 
(actually set of points in state space) de- 
scribed by Equation (9.12) for the case of 
a 23 > 0 and Nj < 0 is examined by the methods 
developed in Chapter 8. In Figure 9.18 , there 
are plots of the three ci.-vs-co. parabolas as 
in Figure 8.1, but with the polarities and 
numbers cf the present equilibrium point. 

The parabolas in Figure 9.18 are degen- 
erate's is alwaystrue whenever an a i3 = 0 
(i = 1, 2, 3) at the equilibrium point. For 
this vehicle, five cases fit this degenerate 
category with only Case I an exception (as 
seen from the table in Figure 9.15) . No real 
stability information is obtained from the 
straight line parabolas other than that (as- 
suming co 2 and co 3 are both positive) a nega- 
tive perturbation of N, and N, drives the sys- 
tem slowly away from equilibrium (as in the 
first 2 parabolas of Figure 8.1), and a posi- 
tive N, andN 3 perturbation drives the system 



-1.0 - 0.75 - 0.5 - 0.25 0.25 0.5 0.75 1 .0 a 23 


(SO') m (- 6 ') (6") (30°) (90°) (sin®) 

Figure 9.17— |<^ 2 Iversus a 23 . 


back (as in the last parabola in Figure 8.1). The important point is that all possible locations 


of equilibrium points have been determined and the computer is used efficiently. For Nimbus, 


the computer responses, starting at many of the equilibrium points, meandered around in almost 



Figure 9.18— Degenerate Stability Parabolas. 
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closed limit cycles, and then converged or diverged. It is apparent that these are marginal 
cases, and second-order effects such as the initial momentum stored in a flywheel (completely 
neglected in this analysis) determine stability. It is emphasized that the only undesirable re- 
sponses that occurred in the hundreds cf thousands examined, were associated with equilibrium 
points! This lends much engineering credence to th epower and usefulness of a thorough equi- 
librium approach. 

It is mentioned in the beginning of this chapter how proposed new computer verified con- 
trol laws for Nimbus II and HI were discarded because of equilibrium points. One of these 
schemes is essentially to raise the yaw torque level N, back up to the same level asN, and If. 
Because of the problem mentioned on page 91 of yaw torque firing mainly due to roll rate, it 
is attempted to cancel out the roll rate information from the gyro compass. This is accom- 
plished by taking the roll horizon scanner signal 4>, putting it through a lead network to get ap- 
proximately 4> + Ktjfc, subtracting the original roll signal 4>, and feeding the result (~ - <£) into the 
yaw channel. That is working with small Euler angles, 

cq % 4> - n o i/i , (9.15) 

where the yaw error term </> is wanted and the roll-rate term <t> is the unwanted coupling. There- 
fore, & is generated so that it can be subtracted from Equation (9.15). At the same time, since 
the gyro-compass signal couples the roll and yaw channels, and since this signal contains five 
times more roll rate than yaw rate, why not feed its output signal into the roll loop? These new 
control laws logically sound good, and they are very good for small angles. As it turned out, 
the computer results via many runs also indicated that the system seemed good, or at least 
better than the original system for large angles. However, the same initial conditions were 
used to compare both systems. This new system is now viewed through an equilibrium approach. 

The simplified control laws are given in Figure 9.19. The steady-state control-law equa- 
tions are 


i - 


£ 3 = V 1 + V 3 


(9.16) 


where v t , v x ' , and v 2 are statically shown in Figures 9.10 to 9.13. v t ' is the same as v x except 
that the saturation level is 30 degrees rather than 12 degrees, 

+ j ^ for |- 9 a j i + 1.9oj 3 J <27.6 deg 

V3 * 2^6sgn(-9co 1 + 1.9^) for |_ 9^ + 1.9a> 3 | >27.6 deg , (9.17) 

(ally’s are in degrees/sec), and the N.'s (i = 1,2, 3) are bang-off-bang functions of the £j as 
shown in Figure 9.8b. 
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Figure 9.1 9-Proposed Control Laws. 

This system has 42 distinct cases which are equilibrium points (Reference 27). The cases 
are established by taking all 27 combinations of the following torque laws: 


N, 


o, 

hi 

< 4 deg 

N. 

= 

0 , 

1 e 2 | 

< 6 deg 

N, 

= 

0 , 

| C 3'| 

< 2 deg 

N ! 

= 

0.35, 


4 deg 

N. 

= 

0.38, 

e 2 > 

6 deg 

N. 

- 

0.19, 


2 deg 

N, 

= 

-0.35, 

< 

-4 deg 

N, 

= 

-0.38, 

C 2 - 

_ 6 deg 

N, 

= 

-0.19, 

e 3 — 

“2 deg (9.18) 

01 

Ni 

<0.35, 

e i 

= 4 deg 

01 

n 2 

10.38, 

€ 2 

- 6 deg 

01 


1 A 

O 

H-* 

LO 

£ 3 

= 2 deg 

-o. 

35 

<Nj <0, 

e i 

= -4 deg 

-o. 

38 

<N 2 10, 

6 2 

- -6 deg 

-o. 

19 

1N 3 10, 

e 3 

- - 2 deg. 


Sometimes graphical techniques are required to solve the algebraic system of relationships 
unless a computer is used for the tedious task. Instead of all of the points, for illustrative 
purposes only one point is investigated. This point is categorized under Case 111 in Reference 
27, and is one of the "holes" in a computer-derived "stable region" mentioned in the beginning 
of this chapter. It is a case for which 

Nj = 0.35 ft-lb £j > 4 deg (full On) 

N, = 0.38 ft-lb e 2 > 6 deg (full on) (9.19) 

0 < N 3 < 0.19 ft-lb e, = 2 deg (duty-Cycling) . 
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Equations (9.7), (9.16), (9.17), and (9.19) are all simultaneously satisfied, and they reduce to: 


w 2 "3 " 

~10. 9(deg/ s) 2 


(9.20a-l) 


+9.52(deg/s) 2 


(9.20a-2) 

oi t o) 2 = 

-130N 3 (deg/s) 2 , 0 < N, 

< 0.19 ft-lb 

(9.20a-3) 


"3 a 23 ~ ^2 a 33 


(9.20b-l) 


03 1 3 3 3 _ ^3 3 1 3 


(9.20b-2) 


"2 3 13 “ “l 3 2 3 


(9.20b-3) 


2 , 2 , 2 _ . 
a 13 + a 23 + a 3 3 _ 1 


(9.20c) 

e l = 

V 1 ( a 1 3 * a 2 3 ) ~ V 3 ("I- " 3 ) 

> 4 deg 

(9.20d-l) 

£ 2 = 

V 2 ( a 1 3 ’ 3 23 ) > 6 de S 


(9.20d-2) 

£ 3 = 

V 1 ( a > r ' a 23 ) + V 3 (" 1 - " 3 ) 

= 2 deg , 

(9.20d-3) 


where v x ', v t and v 2 are given graphically, and v 3 is given in Equation (9.17). 
Equations (9.20a-l, 2) combine with (9.20b-3) to yield 


til - ttl - - - 1.145 , 

a i3 

which plots as a straight line in the a 13 -a 23 plane. 

Equations (9.20a-2, 3) are combined with (9.20b-l) and (9.20c) to yield 



or 


a 


2 

13 


+ 



0 . 00536 \ 

n 3 2 / 


*23 


which is a family of ellipses (0£N 3 < 0.19) in the a 13 -a 23 plane. 


(9.22) 


For the case where N, is full-on positive, negative, or completely off, Equation (9.22) 
represents a single (possible degenerate) ellipse, and has in general eight common solutions 
with Equation (9.21), four each for positive and negative a 33 . Four of these cases are elimin- 
ated by the polarity table, Table 8.2. For the remaining four, Equations (9.20a) are solved ex- 
plicitly for the three cJ s, v 3 is computed, and the control-law inequalities (9.20d) are checked. 
This procedure yields all the equilibrium points. When none of the N's are duty-cycling, the 
procedure is simple and straightforward, since the w's can be determined. It is when one or 
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more of the N's have any value (within bounds) such as this illustrative case, that the pro- 
cedure is more involved. 

The constraint Equation (9.20d-2) permits one to plot the straight line cf Equation (9.21) 
on the graphs of Figures 9.12 and 9.13. The only valid regions are approximately 


0.111 a 13 < 0.66, a33 > ° 

0.04 < a 13 < 0.66, a 33 < 0 , 


(9.23) 


such that e 2 is greater than 6 degrees. Using Equation (9.21) again as a straight line, but this 
time plotted on Figures 9.10 and 9.11, one can graph v 4 and v t ' as a function of a 13 . These are 
shown in Figure 9.20, but only for a 13 positive because of (9.23). The shaded regions are the 
indeterminate regions. No value is defined for a 13 > 0.66, in order not to violate Equation (9.21), 
or (9.23). 

Since the set of Equations (9.20a, b, c) has one more variable than independent relation- 
ship, it is possible to express all the rates and direction cosines in terms of N 3 . Manipulating 
Equations (9.20a) yields 


(9.52) (-130) 
(- 10 . 9 ) 


N„ 


(- 10.9) (-130) 
(9.52) 


N a 


(-10.9) (9.52) 
(-130) N, 


(9.24) 



D(N 3 ) = <-10.9) 2 {9.52) 2 

+ [(-10.9) 2 + (9.52) 2 ] (-130) 2 N 3 2 . (9.25b) Figure 9.20-v, ,v, ' versus a 13 . 
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Hence, once N 3 is known, the equilibrium point is defined. To obtain the values of N 3 , the 
deadband Equation (9.20d-3) is used in a graphical solution. That is, 

V 1 ( a i3> a 23 ( a u)) ~ 2 “ V 3 y'l ( N 3( a 13))’ ^3 ( N 3 ( a 13))^ * (9.26) 

From Equations (9.20a-2) and (9.24), 


V 3 = ~ 9w i fc 


"3 [i q{ - 

h n K)] / ( 


_ (9) (- 130) 
a (-10.9) 3 


( 9 . 27 ) 


But, from Equation (9.25), realizing N, >0, then 


N 


3 



(-10.9) 2 /(-130) 2 

- [l +(-10.9) 2 /(9.52) 2 ] 


(9.28) 


Hence, Equation (9.27) becomes 


v 3 ' V9.52 [l . 9 f (a 13 ) - 9/ f (a 13 )J [sgn , (9.29a) 

where 

1/4 

f(a 13 ) = t l/a i 2 3 - l-(-10.9) 2 /(9.52) 2 ] . (9.29b) 

From Equation (9.17), 

fv 3 ' for jv 3 '| < 27.6 

Vs [^27.6 sgn (v 3 ') for |v 3 '| > 27.6 , (9.30) 

and (2 - v 3 ) is plotted in Figure 9.21 for positive a 13 . Also, v t (a 13 ) from Figure 9.20 is 
plotted in Figure 9.21, and Equation (9.26) indicates that the intersections of the two curves 
are the desired equilibrium points. Thus, there are two determinate-region equilibrium points 
for a 33 > 0, and a range of indeterminate-region points for a 33 <0. All these points satisfy 
Equation (9.23), which ensure satisfaction of (9.20d-2). Equation (9.20d-l) is clearly satisfied 
by Vj' being sufficiently large. Hence, all the equilibrium points corresponding to (9.19) are 
determined. Summarizing the values for the two right-side-up (i.e., positive a 33 ) points, 
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Equations (9.28), (9.24) and (9.25) yield, for 
the first point. 

a 13 = 0.122 


a 23 = -0.139 

(9.313) 

a 33 = 0.983 

(9.31a) 

= -1.08 deg/sec 


“2 ~ 1 -24 deg/sec 

(9.31b) 

“3 ~ -8.79 deg/sec 

(9.31b) 

N, = 0.35 ft- lb 


N, = 0.38 ft-lb 

(9.31c) 

N, = 0.0103 ft-lb 

(9.31c) 

(i.e., N, has a 5.4 percent duty cycle) 




Figure 9.21— v, and 2 - v 3 versus a 13 


Vj - 8.9 deg 
Vj' = 8.9 deg 

v 2 = 7 deg 

Vj - -6.9 deg 


(9.31d) 


e j _ 1 5.8 -deg > 4 deg 

e 2 = 7 deg > 6 deg (9.31e) 

e 3 ~ 2 deg , 

and for the second point, 

a J3 = 0.373 

a 23 = -0.428 (9.32a) 

a 33 = 0.823 
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2 . 08 deg/see 


co 2 - -2.38 deg/sec (9.32b) 

6J 3 ~ 4.58 deg/sec 

N x = 0.35 ft-lb 

N 2 = 0.38 ft-lb (9.32c) 

N, = 0.0380 ft-lb 

(i.e., N 3 has a 20 percent duty cycle) 

Vj = 12 deg 

Vj' - 27.5 deg 

(9. 3 2d) 

v 2 = 12 deg 

v 3 = -10 deg 


e j - 37.5 deg > 4 deg 

e 2 = 12 deg > 6 deg (9.32 e) 

e 3 - 2 deg . 

The first equilibrium point, which is described by Equations (9.31), has a relatively high 
yaw rate, co 3 = -8.8 degrees/sec. The design goal was a three-dimensional rate envelope of 
6 degrees/sec, within which all initial attitudes produced stable responses. Hence, this point 
is not of much concern. However, the second point, which is described by Equation (9.32), is 
definitely of much concern, since the three rates are in the region of interest. Also, the atti- 
tude cf the spacecraft is such that the yaw axis x 3 is only about 35 degrees from the local 
vertical y 3 , that is, a 33 = 0.823. Previous to the calculation of this point, the region in six- 
dimensional state space around the point appeared to be a stable region as a result of many 
computer runs. For instance, initial conditions of 

nq ~ 2 deg/sec 

^ 2 = - 2 . 5 deg/sec (9.33) 

00 3 - 4.5 deg/sec 

were run for 22 strategically located positions on the direction-cosine unity sphere. All the 
responses were well-behaved. The initial rates were varied in ±0.25 degrees/sec step 
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increments such that almost all combinations within the 6 degrees/sec rate envelope were 
tried, each with 22 starting positions, resulting in hundreds cf thousands of computer runs. 

All those near the rates of (9.33) were stable. However, when the initial conditions of (9.32a, b) 
were used (after the equilibrium analysis was performed), the system meandered in an almost 
periodic fashion. The numbers read from the computer were 

a 13 % 0.36 ± 0.07 

a 23 % -0.41 + 0.03 (9.33a) 

a 33 % 0.84 ± 0.02 

co l % 1.95± 0.15 deg/sec 

% -2.4 + 0.2 deg/sec (9.33b) 

&> 3 4.9 ± 0.7 deg/sec 


and 


N, full-on positive 
N 2 full-on positive 

N, on positive with about 25-percent duty cycle. (9.33c) 

This continued for a long time, using up much precious gas; then finally diverged unstably when 
the secondary effects (neglected in any paper analysis) grew sufficiently large. (A typical run 
is Number 181 in Reference 33.) Thus, here is an unstable "hole" in a "stable region." It 
corresponds to an equilibrium point which is not at all obvious, since the sensors have non- 
zero (in fact saturated) outputs, and the gas is firing. As evidenced in this chapter, it is not 
easy to locate these types of points for a highly nonlinear system. 

To conclude this chapter, an advanced Nimbus control system design (Nimbus D whose de- 
sign is just being completed and implementation commenced) is briefly discussed. Emphasis 
is on the equilibrium approach as a design tool. At the outset, the design goal was a system, 
stable in the large, for rates up to 6 degrees/sec in all axes for any attitude. The system was 
to be as simple (hardware-wise) as possible. The previous Nimbus systems all contain only 
one rate gyroscope as a sensor. However, no control scheme could be devised (and some pro- 
posed were quite ingenious) for which no equilibrium points exist in the state-space regions of 
interest. Cti the other hand, it is very easy to design a system with three ideal rate gyros for 
which complete rate stability can be mathematically demonstrated. For example, if one uses 
variable linear thrusters, the well-known control laws are 

N, = -Ki i = 1, 2, 3 , (9.34) 
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where >0, and <a i are the sensed rates. To prove stability, one of many possible suitable 
Liapunov functions that is chosen is the rotational kinetic energy 

v = T J i «i + \ I 2 + J h <4 ’ ( 9 - 35 ) 

which, of course, is positive-definite in “j. The time-derivative is 

V = + I 2 w 2 ^ 2 + 1 3 , (9-36) 


which, when Equations (7.50a) a nd (9.34) are substituted in, is 


V " N, co t + N, a> 2 + N, o> 3 


- “Kj ajj 


2 - K 2 co 2 - K 3 ^ 3 2 


(9.37) 


Equation (9.37) is obviously negative-definite. That is, the control laws cf (9.34) insure that 
rotational kinetic energy is always removed from the rigid body until it comes to rest; but it 
comes to rest at any arbitrary attitude. If one defines ideal position sensors that measure a 13 
and a 23 , then it is possible to evolve control laws which lend themselves to a mathematical 
proof of assymptotic stability in the large, for both rates and positions, via Liapunov functions. 
As an example, from Reference 24, Section 4.3, the following control laws are used: 


N i 


_K l"l - C 2 a 23 


N, - -K 2 a> 2 + C 1 a 13 
N, = -K 3 w 3 . 


(9.38) 


A suitable Liapunov function is 


[li K + a 23 ) 2 + I 2 (o> 2 - A 2 a 13 ) 2 + I 3 a > 2 + B (1 - a 33 ) 2 + (B- I 2 A/) a, 2 , +(B- I, A 2 ) a/J , 


where 


(9.39) 


B = A, K, + C 2 = A 2 K 2 + Cj 

Kj > Aj Ij > 0 

K, > A 2 I 2 > 0 
Cl > 0 
c 2 > ° 

jAj Ij + A 2 (l 3 - I i)J 2 [A 2 I 2 + Aj (l 3 - I 2 )J 

K 3 - 4(K 1 -A 1 I 1 ) + 4(K 2 - A a I 2 ) 


(9.40) 


108 



since V is then positive-definite and V (proof in Reference 24, Section 4.33) is negative-semi- 
definite in co v a 2 , co 3! a 13 , a 23 , and 1 - a 33 . But as demonstrated earlier in this chapter, sen- 
sors are far from ideal and do not yield such things as direction cosines. Also it is more ef- 
ficient to build nonlinear controllers such as on-off thrusters. Hence, a Liapunov or similar 
approach, even with the aid of large-capacity digital computers, has not been fruitful up to this 
time. 


For Nimbus D, it was decided that the approach should be a design that eliminates all (but 
the desired null) equilibrium points from the state-space region of interest using a minimum 
number of reliable components. This was accomplished with a low-cost, non-inertial quality, 
spring-restrained gyro (in addition to the high-accuracy gyro-compass used for yaw position 
information in the fine-accuracy pointing mode). For damping, electronic lead networks are 
employed. The thrusters are driven by pulse-frequency, pulse-width modulators as shown in 
Figure 9.22. These are implemented by placing a passive lag network in a feedback path around 
the usual threshold detector. For a DC input signal x, the output y settles down into a pulse 
train with the width and frequency of the pulses functions of x. The steady-state duty cycle of 
Y, defined as on-time divided by on-time plus off-time, is plotted in Figure 9.22 as a function 
of x. If y is thought of as the driving signal 


to a solenoid for an idealizedpneumatic sys- 
tem, then the graph in Figure 9.22 is one cf 
normalized torque versus x. The advantages 


_FL=JT 

in 'off 


of this control circuit are many. The dead- 
band and full-on characteristics still exist 
asin a pure bang-bang controller (i.e., K = 0) 
for small and large error signals. In addi- 
tion, a linear range exists, since t on +t off 
is small compared to the time constants of 
the overall system loops. That is, the 
steady- state duty cycle acts like a constant 
torque. The loop provides damping, since a 
lag in the feedback path is in effect similar 
to a lead in the feedforward path. (Another 
point of view is that for large T the integral 
of torque, i.e., rate damping information, is 
fed back.) The circuit is analog, simple, re- 
liable, and low-power-consuming. (For a 



flywheel controller, the pulse modulator al- 
lows the use of a very efficiently wound 
2-phase servo motor, with an essentially flat 


ton + toff 
In 


das) 

x-X-KY \ / x-aX V 
x-aX-KY J V x-X / 


torque- speed curve. It is used in an on-off 
fashion, so that the power efficiencyis signi- 


ficantly increased and a linear drive ampli- Fjgure 9.22-Pulse-Frequency, Pulse-Width 

fier is not needed. For analysis purposes, Modulation Control. 
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the system acts as if it is linear for fine control.) The dynamic operating range of the torquer 
is greatly increased. For initial stabilization a large torque is available for quick energy re- 
moval, and yet, for fine unloading of the flywheels, extremely small torques are also available. 
That is, the minimum on-time for a single pulse can be made arbitrarily small, down to the 
limitations cf the hardware, 


minimum pulse width time 


T In 



-a)X 

KY 


(9.41) 



The small pulse time allows the system to 
obtain a final limit cycle that has extremely 


The control laws that are designed for 
Nimbus D are simplified to only essentials 


n. 


GYRO 


o.i°/s 0.257s 

1 ft-lb 


I 3 =145 
SLUG-FT 2 


PULSE MODULATOR 


shown in Figure 9.23. The yaw torque is as 
large as roll and pitch, but is controlled by 
a rate gyro down to 0.1 degree/sec. (The 
gyro-compass still feeds the fine-loop fly- 
wheel which is not shown.) The pulse mod- 
ulators are represented statically. The roll 
and pitch information is displayed graphi- 
tally. However, this time the horizon sen- 
sors and attitude processing is such that no 

electrical clamping is necessary, thus making the error information symmetrical and deter- 
minate, yielding all attitude control. Saturation limits of aZlthe electronics and processing are 
above the full-on levels cf the pulse modulators, and thus do not affect the analysis. 


YAW 


Figure 9.23— Nimbus-D Control i Laws. 


If the assumption is made that the deadband in the yaw loop of Figure 9.23 is zero rather 
than 0.1 degree/sec, then the yaw torque drives oj 3 to zero. This eliminates the cross-coupling 
torque terms in the roll and pitch Euler dynamic Equations (7.50a). Hence, at all equilibrium 
points, Nj and N 2 are zero. Referring again to the control laws in Figure 9.23, and e 2 (and 
thus the roll and pitch computed error signals) lie within the J±5 degree deadbands. For the 
Nimbus D scanners, this occurs only around the desired equilibrium point and at a 13 ±1. 

Hence, aside from the designed null point, the only equilibrium regions that exist are two mirror 
images at approximately: 


a i 3 “ ±x 


= 0 


*23 


*33 


co 2 = ° 

o > 3 = 0 


(9.42) 
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There are no others. These are trivial ones, in that they do not use any gas and the sensors 
have near- zero outputs. They occur when the spacecraft is rotating with an arbitrary rate 
about its roll axis, x x , which is aligned to within about 20 degrees of the local vertical, y 3 . 

Thus one scanner views only the earth and no sky, and the other scanner views only the sky and 
no earth. The earth horizons cannot be distinguished and zero errors are developed. (In this 
attitude the yaw gyro-compass still measures the misalignment between the spacecraft and 
orbit pitch axes, If there is any error, the yaw wheel speeds up and drives the roll spin axis 
off the earth. Thus, (9.42) represents unstable equilibrium points and the system does not long 
remain at the points.) In these attitudes, the solar paddles still fully track the sun, thus not 
depriving the satellite cf any power. If the system gets into these equilibrium conditions, they 
are safe, short lived, and acceptable. 

Much confidence was gained for the Nimbus D control system design by the elimination of 
equilibrium points via driving w 3 to near zero. Many computer runs have been successfully 
made which further solidify one's confidence. However, complete confidence is still out of 
reach, and much work has yet to be done. At present, one must wait for a flight to point out 
any remaining flaws in the attitude control system. Then, using actual telemetered performance 
data, the engineer updates his computer model, and if necessary, betters the system design for 
the next time (as well as figuring out what, if anything, went wrong in orbit). 
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10. SUMMARY AND CONCLUSIONS 


It is shown in Part I, that if a finite bounded limit cycle exists for the 
dimensional system 


general n th - 


*(t> = f(x(t)) , (10.1) 

then an equilibrium point lies inside the limit cycle. That is, 

X im < X ie < X iM • i = 1 , 2, . . . , n (10.2) 

where x e is an equilibrium point satisfying the relationship, 

0 = ?(*.) . (10.3) 

an d' x im 5 x iM are l h e minimum and maximum values of the x(t) = x(t +r) limit cycle solution. 

If no equilibrium points exist, then no limit cycles exist. Equation (10. 2,)represents a physical 
system such that there exists a unique continuous solution x(t) for every initial condition x(0). 

Furthermore, if Equation (10.1) is linear and admits periodic solutions, then one and only 
one equilibrium point exists, and it is at the centroid, 

X e = ' ( 10 - 4 ) 

that is, for periodic solutions x(t) = a(t +r) ,the equilibrium value of each state variable x i is 
the average cf the periodic function (t); 

i = 1,2, ... ,n (10.5) 

In Part II, an equilibrium analysis technique is formulated and illustrated for complex at- 
titude control systems encountered in practice. Since the stability of these systems cannot in 
general be determined by means cf a theoretical analysis, a detailed computer simulation is 
employed. Various combinations of system parameters and initial conditions require a huge 
number of computer runs to gain sufficient confidence in the stability cf the system throughout 
the state-space regions of interest. Using an equilibrium approach significantly raises the con- 
fidence level and lowers the number cf required runs. Troublesome behavior in the physical 
system is linked to limit-cycle behavior in the mathematical formulation, which is shown to be 
associated with equilibrium points. Hence, instead of utilizing auniform or random n-dimensionai 
grid of computer initial conditions, a grid is concentrated around a few calculated points. 


/*t+T 

X. - X. = I : 

16 la Jt 


X (t)dt 
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As discussed in Part II, completely general systems can be handled, but for illustrative 
detail, general two-axis stabilized attitude control systems are considered. They are designed 
to align and hold a spacecraft principle axis to a fixed axis with zero rotational rate (but at an 
arbitrary angle about the axis). The governing equations are 



N. 


+ TT K> "j - “k- a i3> a j3> a k 3 ) 


(10.6a) 


i3 


"k a j 3 


a k3 


(10.6b) 


3 

YL a * 3 = 1 (10.6c) 

m = 1 

where i, j, k are the cyclic indices 1,2, 3. The symbols are defined in Chapter 7; the N. 
torque functions are defined by control laws for a particular attitude control system. Equi- 
librium points are the algebraic solutions cf Equations (10.6) (with all the derivative terms 
zero) that simultaneously satisfy the control laws. The following Equations are derived from 
10 . 6 : 


R 


N i N k 

N. R.. R. . 

i i j ki 


jk 


(10.7) 


or 


Table 10.1 
(I 1 >I 2 >I 3 >0) 



2 

ai3 


N. 


a j3 a k3 ®jk 


(10.8) 


where R £j = l i - . If the moments of inertia 

cf a rigid body are specified, polarity con- 
straints are imposed on the torques and 
direction cosines, so that the right-hand sides 
of Equations (10.7) and (10.8) are positive. 
(These are illustrated by the last six columns 
of Table 8.2 in Chapter 8, and repeated here 
as Table 10.1, where I, > I 2 > I 3 > 0 is as- 
sumed for definiteness.) Equation (10.6b), 
with a i3 = 0 , provides the polarity constraints 
on the rates (with respect to the direction 
cosines). Hence, all equilibrium points sat- 
isfy the polarity relationships given in Table 
10.1. The table is a useful tool to help 
determine equilibrium points. 
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From Equations (10.6), the following equations of motion are derived 


a i3 + 


0 )^+ 0 ), 


2 **})■' 




3 _d 
2 a i3 dt 




^ 2 ) +a j3 (V^j) + < *x“j] _a k3 H ("j (10.9) 


and 


N. 

1 

T7 

1 


. R jk a j 3 a k3 , 

+ as. 

I. a . 2 

1 A 3 









( 10 . 10 ) 


If the state of the system is sufficiently near an equilibrium point such that the angular rate is 
slowly varying with time, then Equation (10.9) approximately reduces to the three uncoupled 
equations , 


a x3 + " 2 ( a i3 “ a x3 equib) “ 0 (10.11) 

where co is the magnitude of the (assumed constant) angular rate vector. From the solution to 
Equation (10.11), it is clear that if a i3 is initially near its equilibrium value, then it tends to 
remain near that value. On the other hand, if the direction cosines are assumed constant near 
an equilibrium point, but the cx>.’s are allowed to vary, then Equation (10.10) reduces to three 
uncoupled equations in 


N. 

I. 

1 


t 


R ik a j3 a k3 


I~ a 


2 

i 3 


( 10 . 12 ) 


Thus, if «. is initially near an equilibrium value, it moves further or closer depending on the 
value of co.. Equation (10.12) is plotted in the 3 phase planes shown in Figure 10.1. The curves 
are parabolas that are stationary with respect to time for constant N i . The intersections cf the 
parabolas with the w. axis (i.e., co. = 0) are the equilibrium points of the system given by 
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Equation (10.8) . The existence of these Points is determined by sgn (NL R jk a j3 a k3 ) , and the 
stability by sgn (N t equib ) = - sgn (R jk a j3 a k3 o> iequib ) as is indicated in Figure 10.1. 

Once the equilibrium points for an attitude control system are computed, the polarities of 
N. a), non-rigorously determine the stability of the points. The stability parabolas shown in 
Figure 10.1 also yield a relative "feel" for the degree cf stability or instability by the slopes 
of the parabolas at the equilibrium points. Also, in Chapter 8, a linear analysis is performed 
to show how stability does not directly depend on the equilibrium values cf the direction cosines. 

Much confidence is gained by a stability analysis of attitude control systems via this equi- 
librium approach as illustrated by some practical applications in Chapter 9. The design of the 
Nimbus D control system that eliminates equilibrium points by control laws which drive a> 3 to 
zero, is an example of how the technique is used as a design tool. However, many computer 
runs are needed, and complete confidence is still beyond reach. Much work remains to be done 
in this area. It is still the actual flight that is the best test for discovering any remaining flaws 
in an attitude control system. Then, using actual telemetered performance data, the engineer 
updates his computer model for the next flight system. 

The attitude control system field is wide-open, interesting, and challenging. Each system 
is unique unto itself, since no general, powerful, analytical technique exists. Future designs 
may be based on nonlinear methods (such as a Fiapunov approach as in Equation (9.39)), and 
these methods may be utilized to build optimum, adaptive, or learning attitude control systems 
via small airborne digital computers. That is, the engineer designs the desired system on paper 
and then solves a synthesis problem, as opposed to the present method cf solving an analysis 
problem for a system designed mainly on the basis of engineering intuition and experience. 
However, in this author's opinion, the synthesis approach is a long way off. The fruitful direc- 
tion for the immediate future is combining semi-analytical tools with the power cf a computer 
simulation. 

Goddard Space Flight Center 
National Aeronautics and Space Administration 
Greenbelt, Maryland, January 4 , 1967 
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Appendix A 

Dynamic Direction Cosine Constraints 


In Chapter 8, from Equations (8.10) and (8.11), that is, 

a i3 = a k a i3 ~ a k3 


(A.l) 


and 


3 

(A. 2) 

i=i 

the complex dynamic relationships of Equation (8.13) are derived. The i, j, k are the cyclic 
indices 1, 2, 3. As in Chapter 7, the angular velocity vector is defined as 



" - Wj Xj + C0 2 X 2 + 0 ) 3 x 3 

w = / ^! 2 + co 2 2 +CO 3 2 , 


(A. 3) 


where x 15 x 2 and x 3 are unit vectors along the body principle axes. y 3 is the unit vector along 
the inertially fixed control reference axis such that, 


a i 3 ~ * y 3 ■ ( A - 4 ) 

When zj is a constant vector with respect to time (i.e., is = o), then Equation (8.13) reduces to 
(8.14), that is. 


a i3 + 0)2 k i3 - 0 ■ (A. 5) 

The solutions to (A. 5) are stated in Equation (8.16), that is, 

a i3 It) = Aj^+Bj^sinwt+C. coswt . (A. 6) 

It is also stated in Chapter 8 that the A i; B. and Cj are constants, but not arbitrary constants. 
That is. Equations (A.l) and (A.2) define only two first-order differential equations. In this 
Appendix, the relationships between the A's, B's, and C's are explored. 

Differentiating (A.6) yields 


coB. cos wt - cdC. sinoj. t 


(A. 7) 
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Equating like terms of (A. 7) and (A.l) yields 


0 

= “k A j 

~ w f A, 


coB. 

A 

= “k c i 

- CO j Q 

(A.8) 

i 

= “k B j 

~ w j B, ■ 



1' the A's, B's, and C’s are considered as nine variables, (A.8) is written in matrix form as. 


0 

“3 

~“2 

0 

0 

0 

1 

1 

0 

0 

0 




s' 

0 

"“3 

0 

“l 

0 

0 

0 

1 

1 

0 

0 

0 


^2 


0 

“2 

-“l 


0 

0 

0 

1 

-t 

1 

0 

0 

0 


A 3 


0 

0 

0 

0 

CO 

0 

0 

0 

—co„ 

3 

“2 


Bi 


0 

0 

0 

0 

0 

CO 

0 

1 

1 

“3 

0 



B 2 

= 

0 

0 

0 

0 

0 

0 

CO 

1 

-1 

1 

■“2 

“1 

0 


b 3 


0 

0 

0 

0 

0 

“3 

■“2 

CO 

0 

0 


Cl 


0 

0 

0 

0 

■“3 

0 

“1 

1 

1 

0 

CO 

0 

C 2 


0 

L° 

0 

0 

“2 

"“1 

0 


0 

0 

CO 


C 3 


0 


The coefficient matrix has a 9 x 9 determinant that is zero for any a>, as expected. Defining 
[P w ] as in Chapter 1, Equation (7.39), Equation (A.9) is rewritten as 



(A.IO) 


where [U] is the 3X3 identity matrix, and (A), (B), (C) are column matrices whose elements 
are A,, B. , C., respectively. Thus (A.8) is expressed by, 


[Pj(A) = 0 
"(B) + [P w ] (C) = 0 
“(C) ~ [P w ] (B) = 0 . 


(A. 11) 
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A good geometric picture is obtained by introducing the vector formalism, 


A = A, x x + A 2 x 2 + A, x 3 

B = B, x l + B 2 x 2 + B 3 x 3 (a12) 

C = qfj+Cjij+tii, , 

as well as s defined in (A .3). Now (A. 11 ) is expressed as. 


A x 

A 

co - 

0 


B x 

CO = 

- c 

(A. 13) 

C x 

CO - 

B , 



where <3 is the unit vector <*>/&>. (w = 0 gives rise to trivial equilibrium solutions, and is not 
considered here.) From (A 13), the following relationships are easily derived: 

A x u - 0 

B • co = 0 

C • CO - 0 

A • B = 0 (A. 14) 

B • C - 0 

C • A = 0 
B = C . 

Therefore, A, B, and Care orthogonal to each other with |b| = j c| , and A is parallel to To, From 
the constraint Equation (A. 2) and the solution Equations (A.6), 


A 2 + B 2 sin 2 cot + C 2 cos 2 o)t t 2(A • B) sin cot + 2(A • C) c os w t + 2(B • C) sinwt coswt 1 . (A. 15) 
Therefore, invoking (A.14), 


A 2 + B 2 sin 2 cot + C 2 cos 2 cot - 1 . 


(A. 16) 
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Thus, the A, B, and C vectors are confined within a unity sphere. That is. 


B 2 = C 2 = 1 - A 2 < 1 , 

A 2 < 1 . 


(A. 17) 


Now, 


A 1 

" ’ = A ’ ?3 = A( A l a 13 +A 2 a 23 +A 3 a 33) 

= x( A j 2 1 A ’ s i n wt * A , C, cos wt +A 2 2 + A, B, sinait + A 2 C, cos wt + A 3 2 +A 3 B, sin wt + A, C 3 cos<^t) 

= -j- (A 2 t (A * B) sin ait + (A * C) cos a>t) 

= A . (A. 18) 

That is, the angle between the A vector (or 5 vector) and the reference vector y 3 is cos' 1 A. 

The geometry is shown in Figure A.l. Since the S3 was assumed constant, the body axes Xj, x 2 , 
x 3 are a set of orthogonal axes fixed with respect to ss anywhere in the unity sphere. These 
axes are rotating at an angular rate w about the <£-axis with respect to y 3 . Thus, a clear geo- 
metrical picture is obtained by keeping S3, A, B. C, S I? x 2 , x 3 stationary in the unity sphere, and 
allowing the plane containing A and y 3 to rotate around A at a rate co. 

To illustrate this dynamic geometry, assume S3 = hkj, i.e., & 2 = o> 3 - 0 and ^ = o>. Thus, 
from Equations (A.8), A, = A, = B, - C, - 0, and C, - - B, and C, = B. . That is. 



Figure A.l — Unity-Sphere Geometry. 


o> - 

A = A, x x 

(A. 19) 

B - B, x 2 + B 3 x 3 

C = - B x 2 + B, x 3 , 

where Aj 2 +B 2 2 +B 3 2 = 1 . Therefore, from 
(A. 6), 

= A, 

a 23 = B, sinwt - B 3 coswl (A. 20) 

a 33 = B, sinwt +B, coswt . 
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formed by y 3 andxj. This plane is rotating 

about the Xj-axis such that, at t - t 0 + -rr/2co, Figure A.2— Dynamics of Unity Sphere, 

x 2 lies in the plane. This example illustrates 

the behavior of a non-accelerating spacecraft rotating about one of its principle axes (x x ) with 
the fixed inertial reference vector at an arbitrary attitude, y 3 . The unity-sphere model is an 
especially useful visual aid for picturing attitude motion when 5 5 is at any arbitrary orientation 
with respect to the x., and the w. are time-varying. 
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"The aeronautical and space activities of the United States sfed&~ be 
conducted so as to contribsae ... to the expansion of husnan kriptvP 
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shall provide for the widest practicable and appropriate dissemination 
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